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ABSTRACT 


As  mathematical  models  become  large  and  complex,  they  become  more 
difficult  to  apply  and  understand.  Repro-model ing  is  a technique  for 
constructing  an  efficient  input/output  approximation  to  a complex  model. 

The  technique  is  founded  upon  an  economical  representation  for  continuous 
piecewise  linear  functions  of  many  input  variables.  Repro-model ing  as 
described  here  is  the  fitting  of  input/output  data  by  a continuous  piece- 
wise  linear  function.  A software  package  called  COMPLIAR  which  accomplishes 
this  is  described,  and  the  theory  underlying  its  use  is  explained.  Instruc- 
tions for  the  use  of  the  software  and  examples  are  included. 


PREFACE 


This  document  describes  an  efficient  computational  procedure  for 
repro-model i ng  a sophisticated  model.  The  procedure  makes  use  of  multi- 
variate linear  regression  theory  to  develop  a continuous,  multivariate, 
piecewise-1 inear  approximation  to  the  relationship  that  exists  between 
the  input  and  output  variables  of  the  original  model.  Specially  designed 
software  (COMPLIAR)  is  described  for  obtaining  the  repro-model  using  repre- 
sentative input/output  data  obtained  from  the  original  model.  Thus,  an 
effort  has  been  made  to  provide  in  one  place:  an  introduction  to  the  repro- 

modeling concept,  an  exposure  to  the  mathematical  principles  involved,  a 

tour  through  a specific  implementation  known  as  COMPLIAR,  and  an  illustra- 

r*  ■ 

tive  example. 

Though  this  document  is  not  a User's  Manual  as  such,  the  reader  wh  • . 

i 

is  primarily  user  oriented  will  find  Sections  4. and  5 most  pertinent,  and 
the  other  sections  may  be  omitted. 
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I.  INTRODUCTION 


1.1  PROBLEM  DESCRIPTION  AND  DESCRIPTION  OF  EFFORT 

The  U.S.  Army  Tank-Automotive  Research  and  Development  Command  (TARADCOM) 
is  concerned  with  improving  the  design  of  wheeled  and  tracked  vehicles  to 
achieve  better  ground  mobility  and  higher  rates  of  survivability.  To  pursue 
these  objectives,  TARADCOM  has  invested  considerable  effort  in  the  develop- 
ment of  a highly  detailed  computer  model  for  vehicle  ground  mobility.  This 
model  was  first  comprehensively  documented  in  [1]  which  describes  the  model's 
status  as  of  July  1971.  This  first  generation  model  is  referred  to  as  the 
AMC  '71  Mobility  Model. 

The  Mobility  Model  has  evolved  as  various  embellishments  were  added  and 
modifications  made  in  order  to  ensure  that  the  model  authentically  accounted 
for  the  principal  factors  affecting  vehicle  mobility.  However,  even  the 

t 1 ■ 

1971  version  of  the  Mobility  Model  is  sufficiently  intricate  to  motivate  the 
use  of  approximations  to  the  model.;  Such  approximations  are  highly  effective 
in  reducing  computer  cost  and  program  size  and  in  simplifying  the  interpre- 
tation of  model  outputs.  This  is  especially  important  when  it  is  desired  to 
use  the  Mobility  Model  as  a design  tool  or  as  part  of  another  model,  e.g.,  a 
survivability  model.  In  these  situations  only  a few  of  the  hundreds  of  param- 
eters which  affect  vehicle  mobility  are  of  immediate  interest.  Yet,  without 
any  model  approximation,  it  is  necessary  to  rerun  the  entire  mobility  model 
in  order  to  investigate  the  effect  of  changing  the  few  parameters  of  interest. 
Extensive  parametric  examination  of  vehicle  mobility  is  impractical  in  this 
situation. 


^TACOM  Technical  Report  No.  11789  (LL143),  The  AMC  *71  Mobility  Model, 
July  1973.  AD766733,  AD766734. 
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To  overcome  this  problem,  TARADCOM  is  exploring  the  application  of 
repro-model i ng  as  a means  of  approximating  the  Mobility  Model. 

Repro-model i ng  [2]  is  an  approach  for  creating  an  efficient  input/output 
approximation  to  a sophisticated  model.  The  repro-model ing  approach  is  based 
upon  obtaining  a simple  mathematical  expression  for  the  desired  input/output 
relationship.  It  is  an  empirical  approach  which  employs  the  original  model 
to  supply  the  requisite  data  (typically  tables  of  model  output  values  for 
a selected  range  of  input  values).  Once  the  Input/output  data  are  at  hand, 
the  idea  is  to  approximate  the  functional  relationship  implicitly  imbedded 
In  the  data  with  a mathematical  expression  which  is  both  economical  to  use 
and  which  reasonably  approximates  the  behavior  of  the  data. 

To  demonstrate  the  potential  utility  of  repro-model ing  in  approximating 
the  Mobility  Model,  TSC,  under  contract  to  TARADCOM,  has  developed  a computer 
program  which  Is  capable  of  repro-model ing  a model  when  the  input  may  contain 
as  many  as  eight  input  variables.  The  computer  program  produces  a multi- 
dimensional functional  approximation  which  relates  the  selected  input  variables 
to  the  output  variable  (the  average  velocity  for  given  terrain  conditions  and 
a given  vehicle  for  the  Mobility  Model).  Once  at  hand,  the  resulting  func- 
tional approximation  may  be  used  as  a design  tool  to  further  study  the  effects 
of  changes  in  the  chosen  input  variables  upon  mobility. 

The  purpose  of  this  document  is  to  describe  the  approach  taken  in 
obtaining  the  repro-model  and  to  completely  describe  the  computer  software. 
Including  instructions  for  its  use. 

^Meisel,  W.  S. , and  D.  C.  Collins,  "Repro-Modeling:  An  Approach 
to  Efficient  Model  Utilization  and  Interpretation,"  IEEE  Transactions 
on  Systems,  Man,  and  Cybernetics;  Vol.  SMC-3,  No.  4,  July  1973. 
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Organizationally  we  first  proceed  to  describe  the  concepts  embedded 
in  the  continuous  Multivariate  Piecewise-LInear  Approximation/Regression 
(COMPLIAR)  program.*  This  is  the  mathematical  approach  which  has  been  utilized 
in  obtaining  the  required  multi-dimensional  approximation.  Next  we  discuss  the 
computer  software  implementation  used  in  COMPLIAR,  including  a step-by-step 
discussion  of  its  use  in  furnishing  a repro-model  of  the  Mobility  Model.  A 
sample  deck  and  the  resultant  computer  output  are  provided  to  illustrate 
the  operation  of  the  software.  (Complete  annotated  listings  of  the  computer 
program  can  be  obtained  from  the  Science  and  Technology  Division  of  the  Tank 
Automotive  Research  and  Development  Laboratory.) 

Together  with  the  actual  computer  program  being  furnished  under  this 

* 

contract,  this  document  should  enable  TARADCOM  personnel  to  exercise  the 
software  and  realize  the  cost  savings  and  interpretative  advantages  of 
repro-model ing.  In  particular,  the  software  may  be  used  to  repro-model  the 
relationship  between  other  input/output  data  sets  involving  up  to  8-dimensional 
input  data  and  a single  output  variable. 

The  software  was  used  on  sample  data  from  the  Mobility  Model  as  a demon- 
stration of  its  use.  One  such  demonstration  utilized  a data  set  consisting 
of  5184  input/output  data  points  describing  the  mobility  of  an  M60  Tank  as 
a function  of  8 input  variables  (4  input  variables  described  terrain  proper- 
ties and  the  other  4 described  vehicle  design  parameters).  The  repro-model 
for  this  test  case  accounted  for  more  than  94%  of  the  functional  variation 
contained  in  the  data  and  required  approximately  60  CPU  seconds  of  execution 
on  a CDC  6400  to  completely  specify  the  repro-model.  The  repro-model,  itself. 

While  one  should,  in  general,  distinguish  an  algorithm  (e.g.,  continuous 
piecewise  linear  regression)  from  a specific  software  implementation  of  that 
algorithm  (e.g.,  the  COMPLIAR  software  package),  we  will,  for  the  sake  of 
brevity,  not  be  precise  in  making  this  distinction. 
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can  easily  be  implemented  on  a hand  calculator  (See  Section  5.2  page  110) 
to  furnish  small  volumes  of  input/output  data. 

1.2  MOTIVATIONS  FOR  REPRO-MODELING 

It  is  worthwhile  to  review  here  some  of  the  salient  factors  which 
motivate  the  use  of  repro-models. 

Complicated  and  intricate  models  that  require  excessive  running  time 
and/or  computer  storage  may  be  replaced  by  far  more  efficient  "repro-models" 
that  are  largely  equivalent  for  certain  purposes  to  the  original  models. 
Repro-models  are,  effectively,  simpler  models  of  such  complex  models.  Repro- 
models  make  practicable  many  model  applications  that,  heretofore,  would  not 
have  been  considered  because  of  exceedingly  large  requirements  on  computa- 
tional resources. 

Typically,  a great  deal  of  thought  and  effort  is  expended  in  the  devel- 
opment of  a complex  model.  However,  when  the  model  is  to  be  used  operation- 
ally, it  often  has  to  be  greatly  simplified  or  used  sparingly  because  of 
excessive  requirements  for  computer  time  or  storage.  The  net  result  is  that 
much  that  the  model  has  to  offer  is  sacrificed— and  the  cost/benefit  ratio 
of  model  usage  appears  to  be  unattractive. 

Repro-model i ng  is  often  an  alternative  to  this  situation.  The  result 
is  that  the  full  power  of  the  sophisticated  model  may  be  brought  to  bear  on 
analysis  of  complex  problems,  and  realistic  solutions  may  be  obtained.  For 
example,  complex  systems  analyses  may  now  incorporate  fully  realistic  models, 
yet  require  minimal  computer  resources. 
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In  the  case  of  field  deployment  of  sophisticated  algorithms  for  signal 
or  data  processing,  the  use  of  repro-models  may  allow  these  algorithms  to 
be  employed  in  real-time  and  within  very  limited  storage  constraints  (e.g., 
within  microprocessors)  without  sacrificing  performance. 

Furthermore,  analysis  of  the  repro-model  input/output  data  yields  an 
understanding  of  what  the  model  implies,  and  consequently  provides  a means 
for  model  verification  or  improvement. 

For  the  Army  Mobility  Model,  the  major  objective  of  repro-model ing 

• $ " 

is  to  accomplish  a reduction  of  model  running  time  and  program  size.  This 
is  a typical  objective  of  repro-model ing.  For  example,  a photochemical  air 
pollution  repro-model  developed  by  TSC  for  the  Environmental  Protection 
Agency  [3]  was  so  simple  that  the  repro-model  output  could  be  Calculated 
on  a hand  calculator  in  a minute  or  so;  while  the  original  model  took  one- 
half  hour,  using  a large  general  purpose  computer,  to  produce  the  same 
results.  Also,  the  simple  form  of  the  photochemical  repro-model' greatly 
facilitated  an  understanding  of  the  degree  of  reduction  in  pollutant 
emissions  required  in  order  to  achieve  a given  level  of  air  quality,  as 
implied  by  the  original  complex  model.  It  also  revealed  weaknessess  in 
the  original  complex  model,  particularly  a large  dependence  of  the  output 
on  assumed  boundary  conditions. 

It  is  often  not  obvious  that  repro-model ing  is  applicable  in  a given 
application  without  detailed  examination.  For  example,  some  less  obvious 
uses  of  repro-model ing  analyses  include  the  following: 

^Horowitz,  Alan,  W.  S.  Meisel,  and  D.  C.  Collins,  "The  Application 
of  Repro-Modeling  to  the  Analysis  of  a Photochemical  Air  Pollution  Model," 
EPA  Report  No.  EPA-650/4-74-001 , December  1973. 
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- Determination  of  which  variables  most  affect  the  output:  the  degree 

of  error  incurred  in  ignoring  (or  not  having  measurements  of)  the 
less  important  variables. 

- Avoiding  on-line  optimization:  a model  involving  an  optimization 

process  can  be  run  off-line,  and  the  relationship  of  input  variables 
to  optimum  values  of  parameters  can  be  repro-modeled. 

- Repro-model i ng  differential  equations:  a time  series  of  values 
generated  by  a model  can  be  analyzed  to  obtain  a simple  difference 
equation,  perhaps  with  a larger  time  step,  which  closely  approximates 
the  behavior  of  the  original  model. 

It  is  hoped  that  the  above  discussion  will  stimulate  the  imagination  of 
those  who  will  have  an  opportunity  to  apply  the  software  described  in  this 
document  to  repro-model i ng  the  Mobility  Model  and  other  complex  models  used 
by  the  Army. 
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II.  CONTINUOUS  MULTIVARIATE  PIECEWISE-LINEAR 

approximation/regression  (compliarT 

As  stated  in  the  Introduction,  the  objective  of  repro-modeling  is  to 
furnish  an  approximate  functional  relationship  between  the  inputs  and  outputs 
of  a complex  system  so  that  the  repro-model  is  more  efficient  to  run  and  eas- 
ier to  interpret  than  the  original  model. 

One  means  of  achieving  this  objective,  insofar  as  the  Mobility  Model 
is  concerned,  is  to  make  use  of  a technique  which  we  will  refer  to  as  Contin- 
uous Multivariate  Piecewise-Linear  Approximation/Regression  (COMPLIAR). 

2.1  COMPLIAR  DEFINED 

Mathematically,  COMPLIAR  is  a technique  for  creating  a functional  form 
which  approximates  the  multivariate  relationship  between  a set  of  input 
variables  and  one  or  more  output  variables,  . 

Thus,  consider  the  system  illustrated  in  Figure  1.  For  the  case  at 
hand  the  "system"  constitutes  the  Mobility  Model.  The  inputs  consist  of  a 
set  of  variables  such  as  vehicle  design  parameters  and  terrain  descriptive 
parameters,  and  the  output  consists  of  a measure  of  the  ground  mobility  of 
the  vehicle  such  as  the  maximum  velocity  attainable  under  a given  set  of 
terrain  conditions. 

It  is  convenient  to  designate  the  input  variables  x, , Xg,...,  xn 
utilizing  vector  notation;  so  that 

x = (x1  — xk — xn)  (2-1) 

describes  the  vector  input  to  the  system.  In  the  present  application  in 
which  the  output  is  the  single  quantity,  speed,  it  is  convenient  to  desig- 
nate the  generic  output  variable  with  the  scalar  y. 
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Figure  1.  Vhe  Army  Mobility  Model 
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Underlying  our  entire  procedure  is  the  fact  that  there  exists  a rela- 
tionship between  y and  x which  we  may  summarize  by  writing 

y = f (x)  (2-2) 

That  is,  there  exists  a rule  for  mapping  an  n-dimensional  point  in  space 
determined  by  the  components  of  £ into  a point  (one-dimensional  space) 
determined  by  the  value  of  y.  The  objective  of  repro-modeling  may  be 
mathematically  stated  as  arriving  at  an  efficiently  executed,  and  accurate 
approximation  to  the  function  f(*)  using  sets  of  values  of  x and  the  con- 
comitant variable,  y. 

i.  L 

We  shall  denote  the  vector  data  point 

/ « (x\  /)  ; l = 1,2,. ..,M  (2-3) 

and  the  set  ■ .j  . - 

jr={z1“*zMJ  (2-4) 

will  be  used  to  designate  the  set  of  M observations  upon  which  the  repro- 
model is  to  be  based. 

Now,  there  are  many  mathematical  techniques  for  approximating  f(*) 
given  the  data^T  For  example,  one  can,  in  principle,  use  an  n-dimensional 
Taylor's  series  which  would  lead  to  a polynomial  approximation.  The  tech- 
nique of  COMPLIAR,  which  we  are  about  to  elaborate  upon,  has  certain  features 
which  render  it  more  attractive  to  the  problem  at  hand  than,  for  example, 
polynomial  approximations.  Among  its  advantages  are:  ease  of  interpretation. 
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better  extrapolation  behavior,  and  a greater  degree  of  parsimony  than  poly- 
nomial expansions  for  approximating  functions  with  regions  of  sharp  curvature. 

Before  embarking  upon  the  mathematical  details  of  COMPLIAR,  let  us 
remark  on  the  properties  of  the  technique  which  are  suggested  by  its  name 
alone. 

The  continuous  aspect  of  COMPLIAR  refers  to  the  fact  that  the  sought 
for  approximation  will  be  a function  which  is  continuous  in  all  of  its  vari- 
ables. This  is  not  always  a requirement  of  repro-modeling.  Indeed,  there 
are  systems  which  exhibit  discontinuous  behavior,  at  least  macroscopically. 
However,  for  the  application  of  repro-modeling  to  the  Mobility  Model,  espe- 
cially with  design  study  applications  in  mind  in  which  one  typically  wishes 
to  examine  effects  of  relatively  small  changes  in  the  input  variables*,  the 
continuity  requirement  seems  eminently  reasonable. 

The  piecewise- linear  aspect  of  COMPLIAR  describes  the  underlying  func- 
tional form  which  is  used  in  order  to  compose  the  final  functional  approx- 
imation. 

A piecewise- linear  function  is  one  which  is  linear  in  subregions: 

1allxl  + a12x2  + - + alnxn  + a1,n+1  for  e X1 

: (2-5) 

aRlxl  + aR2x2  + "•  + aRnxn  + aR,n+l  for  * £ XR 

where  Xj»  Xg X^  are  disjoint  subregions  covering  the  region  of  interest 

* 

for  the  input  variables.  Each  line  of  Equation  (2-5)  is  specified  by 

*x  e X,  Is  the  conventional  mathematical  shorthand  for  "the  vector  £ 
is  in  the  set  Xj." 


n 

(n  + 1)  constants:  the  n coefficients  of  x^,  Xg*..^  plus  a constant  term. 

Once  these  constants  are  specified,  each  line  of  Equation  (2-5)  defines  an 

* 

n-dimensional  hyperplane  in  the  (n  + 1 )-dimensional  space  spanned  by  the 
data  points,  z_  . For  this  reason  we  refer  to  the  constants  in  Equation  (2-5) 
as  hyperplane  coefficients.  The  determination  of  these  coefficients  is 
commonly  referred  to  as  a regression  problem. 

The  idea  behind  COMPLIAR  is  to  suitably  obtain  the  disjoint  subregions 
and  hyperplane  coefficients  so  that  a satisfactory  approximation  to  the  data 
can  be  obtained  within  the  framework  of  the  piecewise-1 inear  component 
structure.  An  example  of  such  a fit  in  the  one-dimensional  case  (y  = f(x-|)) 
is  illustrated  in  Figure  2.  Figure  3 is  an  example  of  a fit  involving  two 
input  variables  (n  = 2). 

The  remainder  of  this  section  is  devoted  to  the  discussion  of  an  effi- 
cient computational  algorithm  for  implementing  COMPLIAR.  The  algorithm 
yields  a continuous  multivariate  piecewise-1 inear  function  over  general 
(implicitly  determined)  subregions  using  a strategy  which  minimizes  the  mean 
squared  error  in  the  final  approximation. 

2.2  IMPLEMENTING  COMPLIAR 
2.2.1  Basic  Formulas 

XL 

Let  x^  denote  the  k— component  of  the  vector,  x_,.  of  input  variables. 

The  vector  £ defines  a point  in  n-dimensional  space  for  each  prescribed 
value  of  the  input  variables.  For  notational  convenience,  it  will  prove 
to  be  convenient  to  augment  the  x vector  with  an  additional  element  of  unity. 

* 

In  the  two-dimensional  case  (n  = 2)  a hyperplane  reduces  to  the  con- 
ventional notion  of  a plane.  In  one  dimension  a straight  line  results. 
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Thus,  let 


X*  = [x1  •••  xk  •••  xn  !] 


(2-6) 


Now  consider  the  vector  of  (n  + 1)  parameters 


V Caqi  V V •"  Vn+i)3 


(2-7) 


Then  the  function 


Hq(i)  - a,-  x'T  - t aqk  xk  ♦ aq>(n+n 


(2-8) 


defines  an  n-dimensional  hyperplane,  and  a comprises  the  vector  of  hyperplane 
coefficients. 

We  wish  to  compose  a function  from  hyperplanes  of  the  form  given  by 
Equation  (2-8),  in  such  a way  that  the  composed  function  is  continuous  in  x. 
Thus,  suppose  we  set  up  a collection  of  Q hyperplanes  via  the  set  of  coef- 
ficients {a_.;  q = 1,  2,...,  Q).  These  coefficients  may  be  conveniently 
— M 

arranged  into  the  Q x (n  + 1)  matrix  shown  below. 


|4-(n  + lH 


I 

A = Q 
"1 


S-1 


Sq 


(2-9) 


Superscript  T above  a vector  or  matrix  will  be  used  to  denote  its 
transpose. 
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Then,  define  the  function 

P(x)  = Max  -CH  (>l) > = P(x;A)  (2-10) 

” l<q<Q  q 

The  function  P(x)  is  called  a P-funotion.  It  has  the  desired  continuity 
property  and  shall  serve  as  a fundamental  building  block  for  arriving  at 
the  final  functional  approximation.  The  fact  that  the  P-function  depends 
implicitly  upon  the  entire  set  of  hyperplane  coefficients  has  been  empha- 
sized via  the  notation  P(x;A)  on  the  far  right  of  Equation  (2-10).  Ulti- 
mately we  will  describe  a procedure  for  adjusting  A to  accomplish  the 
desired  functional  approximation. 

Notice  in  Equation  (2-10)  that  the  maximization  operation  causes  one 
to  "switch"  from  one  hyperplane  to  another  in  a continuous  fashion,  depending 
upon  which  hyperplane  produces  the  largest  value  of  H (x)  for  any  poin*  x» 

M 

In  between  these  switches  one  remains  on  a given  hyperplane.  The  ability  to 
switch  from  one  hyperplane  to  another  permits  us  to  produce  basic  functional 
building  blocks  with  P(_x)  that  are  highly  nonlinear  and  continuous;  even 
though  the  constituents  are  simple  linear  functions  of  x.  An  example  of  a * 
P-function  is  illustrated  in  Figure  4.  This  is  a one-dimensional  example 
composed  of  3 hyperplanes.  The  hyperplane  coefficients  are  specified  by 

* [-.5  -1.5] 

% - to  -1] 

2.3 -[1  -Ij 
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P(x) 


Note:  Solid  curve  depicts  the  P- function 

which  results  from  applying  Equation  (2-10) 
to  the  3 lines  illustrated  above.  The 
dashed  extensions  show  the  continuation 
of  the  component  line  segments. 


Figure  4.  A P-function  in  One  Dimension 


- and  the  P-f unction  is  determined  from 

P(x)  = Max  j[-.5x  -1.5]  ; -l;[x-l]J 

However,  there  is  one  important  limitation  to  the  form  of  P(xJ  and 

* 

that  is  that  it  always  produces  a convex  function.  (This  is  a consequence 
of  the  maximization  operation,  which  happens  to  be  an  easy  computational 
method  of  achieving  the  continuity.) 

Not  all  functions  that  we  shall  choose  to  approximate  are  convex.  We 
therefore  wish  to  modify  our  approach  so  that  we  may  approximate  an  arbitrary 
continuous  function  while  still  preserving  the  simplicity  of  the  functional 
building  blocks  described  above. 

We  may  remove  the  convexity  limitation  by  considering  the  following 

<*  • 

generic  form  for  the  approximating  function: 

• * 

Np 

F(x)  =^w.  P.(x)  + wN +1  (2-11) 

i=l  P 

Thus,  we  simply  weight  the  sum  of  Np  P-functions  using  the  weights 
{w^;  i = 1,...,  Np+1}.  The  use  of  both  positive  and  negative  weights  over- 
comes the  convexity  limitation. 

Our  problem  is  now  to  take  the  data,^T  and  fit  it  with  F(xJ.  To  per- 
form the  fit  we  select  a value  for  Np,  the  number  of  P-functions  which  our 
fit  will  consist  of.  We  also  specify  the  number  of  hyperplanes  which  will 
comprise  each  P-f unction.  Extending  the  notation  employed  in  Equation  (2-9), 

A convex  function  is  "bowl -shaped". 
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we  let  Q.  denote  the  number  of  hyperplanes  in  the  i—  P-function.  The  hyper- 
plane coefficients  for  the  i—  P-function  are  summarized  by  the  matrix  A^, 

/.  \ 

which  has  dimension  Q..  x (n  + 1).  The  rows  of  A.  are  denoted  by  {a^  ; q = 1, 
2,...  Q^>.  The  vector  comprises  the  coefficients  of  the  q^-  hyperplane 
in  the  1—  P-function.  We  then  have,  for  the  i—  P-function  (using  a gener- 
alized form  of  Equation  (2-10)): 


P<(x)  = 


(D  % 

Max  {H  (x)> 
i<q<Qi  q 


= p(x;A.) 


(2-12) 


where  the  notation  P(x;A.)  on  the  far  right-hand  side  of  (2-12)  stresses  the 
dependence  of  P.(xJ  upon  the  matrix  of  hyperplane  coefficients.  A.. 

The  generic' form  of  our  approximation,  which  is  given  by  Equation  (2-11), 

is  seen  to  be  implicitly  dependent  upon  the  weights  {w^;  i = 1 Np+1}  and 

the  hyperplane  coefficients  (A^;  1 = l,...,Np}.  We  shall  refer  to  these 
aspects  of  the  fit,  for  economy  of  notation,  using  the  (Np+l)-dimensional 
weight  vector 

(2-13) 

and  the  overall  coefficient  matrix 


W = w1  w2  * 
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Once  the  values  -of  Np  and  {Q^  = i = l,2,...Np}  are  specified,  the  function 
F(x)  is  specified  fully  by  choosing  A and  W.  To  emphasize  this  dependence 
we  often  write  A,  W).  The  objective  is  now  to  determine  A and  W so 
that  F (?c_;  A,  W)  provides!  a good  approximation  to  the  input/output  data. 
2.2.2  Mean-Squared  Error  Optimization  Criterion 

In  order  to  finally  determine  A and  Wit  is  necessary  to  adopt  some 
criterion  of  goodness  for  the  resulting  functional  approximation.  We  shall 
use  the  optimization  criterion  of  minimizing  the  mean-squared-errcr  (MSE) 
of  the  approximation  for  the  specified  values  of  Np  (number  of  P-functions) 
and  {Q.. ; i = l,2,...Np>  (number  of  hyperplanes  per  P-function). 

For  M data  samples  the  MSE  may  be  expressed  as 


E ■ 


1 

M 


- p(2L£ ;A,w)] 

£*■  1 


(2-15) 


The  minimization  of  E must  be  done  numerically,  and  we  shall  see  that 
the  requisite  computations  can  become  quite  costly  for  large  values  of  Np 
and  Q..  For  this  reason,  we  typically  choose  the  number  of  P-functions  and 
number  of  hyperplanes  per  P-function  to  be  reasonably  small  during  our  first 
fitting  attempts.  Of  course,  if  Np  = 1 we  lose  the  ability  to  fit  non-convex 
data;  this  can  be  a useful  tool  for  discovering  if  the  data  is  intrinsically 
convex.  After  examining  the  quality  and  characteristics  of  the  fit,  it  is 
easier  to  make  judgments  about  further  increases  in  the  complexity  of  the 
fit  to  be  attempted.  Perhaps  surprisingly,  even  Np  =22  or  3 with  comparable 
values  of  Q..  often  produces  a totally  acceptable  fit  because  of  the  inherent 
adjustability  of  the  approximating  function. 
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Insofar  as  the  minimization  of  E Is  concerned  it  turns  out  (see  Section 

* 

2.2.3)  that  for  any  specified  A there  is  a closed-form  analytic  solution 
for  the  optimum  set  of  weights,  W.  The  reason  that  an  analytic  solution 
for  W is  obtainable  Is  that  once  A Is  specified,  the  P- functions  are  com- 
pletely determined.  F(x)  then  depends  only  upon  W,  and  the  dependence  is 
linear.  By  way  of  contrast,  F(x_)  depends  nonli nearly  upon  A,  and  there  is 
no  closed-form  analytic  result  for  the  A which  minimizes  E. 

The  fact  that  W Is  completely  determined  by  A Insofar  as  minimizing  E 
Is  concerned,  implies  that  the  total  number  of  "free  parameters"  that  actually 
benefit  the  fit  is  equal  to  the  total  number  of  elements  in  the  overall  hyper- 
plane coefficient  matrix.  A,  in  Equation  (2-14).  Thus,  the  total  number  of 
free  parameters  in  the  overall  fit  is  equal  to  • 

"p 

Q,  (2-16) 

1=1 


For  example,  suppose  the  problem  at  hand  requires  fitting  the  relationship 
between  4 input  variables  and  one  output  (dependent)  variable.  In  this  case 
n = 4.  Even  the  simplest  non-convex  fit  using  two  P-f unctions  (Np  = 2)  and 
two  hyperplanes  per  P-function  (Q-j  = Qg  = 2)  produces  Np  = 20.  This  fit 
Involves  a numerical  optimization  over  a 20-dimensional  parameter  space  in 
order  to  minimize  E.  In  subsequent  sections  we  shall  describe  a practical 
computational  procedure  for  accomplishing  the  requisite  optimization.  First, 
it  is  in  order  to  set  up  the  optimization  equations. 

* 

As  opposed  to  requiring  numerical  minimization. 
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The  computational  problem  before  us  is  to  minimize  the  MSE,  E,  by  suit- 
ably adjusting  the  collection  of  hyperplane  coefficients.  A,  and  weights,  W. 

Formally  speaking,  we  may  regard  E as  a function  of  A and  W and  denote 
this  via  the  notation  E(A,W).  It  can  be  shown  that  E is  differentiable  with 
respect  to  the  parameters  A and  W.  Thus,  a necessary  condition  for  a minimum 
of  E is  that  the  total  derivative  of  E with  respect  to  all  of  the  parameters 
of  A and  W be  equal  to  zero. 

To  elaborate,  it  is  helpful  to  reassemble  the  elements  of  A and  W into 
one  overall  row  vector  and  to  utilize  the  vector  derivative  or  gradient  con- 
cept. A has  Np  elements  and  W has  (Np  + 1)  elements.  Thus,  suppose  we 
define  the  vector 

|-«-Np-*J-«-Np+l-*-| 

I'  t “ i M ] - Cy,  Y2  "•  Yy  — Y(„f+N|)+1)]  (2-17) 

where  a is  an  Np-dimensional  row  vector  composed  of  the  elements  of  A and  W 
is  as  defined  in  (2-13).  In  other  words,  there  is  a one-to-one  correspondence 
between  the  elements  of  the  row  vector  a and  the  elements  of  the  matrix  of 
hyperplane  coeffi cents,  A.  The  manner  in  which  the  elements  of  A^ are  rear- 
ranged into  the  vector  a is  arbitrary.  Later  we  will  offer  one  possible 
arrangement. 

Now,  strictly  speaking  E may  be  regarded  as  a function  of  the  composite 
parameter  vector  y.  Each  specification  of  the  elements  of  y defines  a point 
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In  the  overall  parameter  space,  r,  which  has  dimension  Np  + Np  + 1.  A mini- 
mum of  E occurs  only  if  the  magnitude  of  the  gradient  of  E (also  called  the 
“norm"  of  E)  calculated  in  the  r space  is  zero.  We  denote  the  gradient  of 

E via  the  notation  . This  is  a vector  defined  as  follows: 

dy 


Vi  V1+“*+ 


9y(Nf+Np+1)  r(Np+Np+1 ) 


3E 

W 


where  r.  denotes  a unit  vector  in  the  i—  direction.  As  indicated  in  (2-18) 

the  first  Nr  components  of  comprise  the  gradient  of  E with  respect  to  the 
r ay  + 

— *\r 

components  of  A,  and  the  last  (ND  + 1)  components  of  comprise  the  gradient 

— gr 

of  E with  respect  to  the  elements  of  W.  In  other  words,  ^ is  calculated 

8E  — 

holding  A wnstant,  and  is  calculated  holding  W constant.  Therefore,  in 
order  for|||jto  be  zero  for  arbitrary  adjustments  in  A and  If  we  require 
separately  that 
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and 


(2-19a) 


(2-19b) 


Using  Equations  (2-11)  and  (2-15)  it  is  possible  to  explicitly  solve 
(2-1 9b)  above  for  the  optimum  W.  In  particular,  let 


9^-j  = ^ -j (ii.  », iLj ) * ^ = l»»»«»Np 


(2- 20a) 


(2- 20b) 


That  is,  6 is  an  M x (Np  + 1)  dimensional  matrix  composed  of  the  P-functions 
evaluated  at  the  various  input  data  points,  to  which  a column  of  l's  is 
adjoined.  From  (2-11)  and  (2-20a  & 20b)  we  can  define 
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F£  = F£(xA)  = ^ wi  P1(xA)  + wN  +1  for  £ = 1 M (2-21a) 

i=l  P 


F=£F1  F2-F£-  V 


from  which  it  follows  that  £ can  be  expressed  in  the  form 


F = W 6' 


(2-22) 


and  E can  be  expressed  in  the  form 


WGT]  jy- WGT] 


(2-23) 


Formally  calculating  yields: 


w • I QL£Tfi-iG]  *§CE.-ri  £ 


(2-24) 


gr 

Setting  = 0 results  in  the  solution: 


w * £ £ [iT  g]  ' * w(A) 


(2-25) 
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where  we  have  emphasized  the  fact  that  the  solution  is  dependent  upon  A 
(through  (3)  via  the  notation  W(A).  It  can  also  be  verified  that  — 5- > 0; 

awr 

so  the  solution  in  (2-25)  does,  indeed,  provide  a minimum  for  E. 

T 

The  solution  for  W given  in  (2-25)  exists  as  long  as  the  matrix  [£  Gj 
is  invertible.  One  requirement  for  this  to  be  true  is  that  the  number  of 
sample  points,  M,  used  in  the  fit  must  be  at  least  equal  to  (Np  + 1).  That 
is,  (2-25)  has  a valid  solution  only  if* 

I 

M £ (Np  + 1)  (2-26) 

Since  the  number  of  P-functions  comprising  the  fit  is  typically  rather  small, 
(2-26)  is  usually  satisfied  in  practice. 

If  (2-26)  is  satisfied  but  [G1  Gj  does  not  have  an  inverse,  it  is 
because  two  or  more  rows  or  columns  of  [G^  G]  are  linearly  dependent.  If 
this  occurs  the  problem  can  generally  be  overcome  by  reducing  the  number  of 
P-functions  and/or  hyperplanes  which  are  being  tried  in  the  fit. 

If  now  remains  to  calculate  r—  . The  magnitude  of  the  y—  component  of 
this  Np-dimensional  vector  is  obtained  using  (2-15): 


2 

M 


M 

5^[y£  - F(xA;A,W)] 

A=1 


9F  (x*') 
3a.. 


U=l,2 NF  (2-27) 


* 

Actually  a more  important  practical  requirement  is  that  M»Np.  This 
ensures  that  the  data  is  not  overfit  by  simply  having  an  unjustifiably  large 
number  of  parameters  in  the  fitting  function.  Note  from  (2-16)  that  M»N,- 
implies  that  (2-26)  is  satisfied.  h 
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In  order  to  proceed  further  it  is  necessary  to  specify  the  correspondence 
between  the  elements  of  a and  the  entries  of  the  hyperplane  coefficient 
matrix,  A. 

Referring  back  to  Equations  (2-9)  and  (2-14)  the  following  correspon- 
dence is  suggested. 

Take  the  elements  of  the  matrix  A.,  (these  are  the  hyperplane  coeffi- 
cients associated  with  the  i—  P-function)  and  assemble  them  into  a row 
vector,  a...  To  do  this  let  a.,  be  assembled  as  shown  below 


1—  row  of 

*1 


2—  row  of 
-i 


q — row  of  last  row  of 


(2-28) 


In  other  words  a^  are  the  hyperplane  coefficients  of  the  q^-  hyperplane 
of  the  i—  P-function  assembled  into  an  (n  + l)-dimensional  row  vector. 

Then  a.  collects  all  of  the  coefficients  of  the  i—  P-function  into  a single 
row  vector  of  dimension  [Q.  • (n  + 1)]  by  adjoining  the  partitions  as  indicated 
in  (2-28)  above. 

Finally  we  form  the  overall  a vector  by  letting 


“'tajsz  ! •••  j 2ii 


OLr 


•a 


•aNf] 


(2-29) 


J.L 

The  vector  a is  therefore  partitioned  into  Np  subvectors  in  which  the  i— 

i.  L 

partition  comprises  all  of  the  coefficients  of  the  i—  P-function.  Equation 

f*  h 

(2-29)  establishes  the  correspondence  between  the.  ]i—  component  of  a and  the 
appropriate  P-function  and  hyperplane  coefficient.  It  is  now  possible  to 
evaluate  (2-27). 


th  • j?. 

Consider  the  i— 1 P-function.  For  any  given  input  data  point  x and  a 

trial  solution  for  we  can  evaluate  (using  the  generalized  form  of  Equation 

(2-8)) 


H^tx*)  = • x,T  ; q = 1,2, ...Q. 


(2-30) 


/■M  o 

We  can  therefore  identify  the  value,  q , for  which  Hv  '(x_  ) is  maximum  over 

u M 

**  fh 

the  range  of  q.  From  the  definition  of  the  i—  P-function  (Equation  (2-12)) 
and  from  the  generic  form  of  F(x)  given  in  (2-11)  it  follows  that  aFfx^J/aa^ 
(which  appears  in  (2-27))  is  identically  zero  unless  a is  one  of  the  elements 


of  the  vector  a_ 

«(1); 

K: 


Moreover,  from  (2-30)  we  obtain  for  the  k— element  of 


SF(x£) 


3a 


rrr 

qok 


= wi  xk 


k = 1,2, ...n 


(2-31) 


For  other  values  of  q f qQ,  the  derivative  is  zero. 

Equation  (2-31)  holds  for  all  values  of  i.  Thus,  to  evaluate  the 
quantities  3F(x_A)/da  appearing  in  (2-27)  one  proceeds  as  follows: 
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o 

(1)  Choose  a data  point  x. 

(2)  Calculate  (x*)  for  i = 1 Np 

(3)  For  each  i above,  identify  the  value  of  qQ  for  which 

H*,1’  (**)  > H<*>  (x*)  ; q = 1,2 Q, 

This  identifies  which  row  of  A.  will  contribute  a 
nonzero  result  to  the  final  derivative.  More  specifi- 
cally it  identifies  which  partition  of  a.  In  (2-28) 
will  contribute  a nonzero  result. 

(4)  Utilize  the  association  between  and  the  partitions 

set  forth  In  (2-28)  and  (2-29).  This  identifies 
0 

the  components  BF(x  )/aa  which  are  identically  zero. 

r 

In  addition.  Equation  (2-31)  provides  the  values  of  the 
nonzero  entries. 


To  obtain  the  desired  gradient,  we  must  evaluate  (2-27)  by  repeating  the 

o 

above  steps  (1)  - (4)  for  each  of  the  M data  points  {x  ; SL  - 1,...,M}  . 

After^havlng  evaluated  the  components  of  using  (2-27)  and  using  (2-24) 

3E  ~~ 

for  -gg  , the  norm  of  the  gradient  can  be  calculated  in  the  conventional 


manner: 


Observe  in  the  above  discussion  that  the  resulting  gradient  is  func- 
tionally dependent  upon  W and  A.  However,  we  have  already  shown  (Equation 
(2-25) ) how  to  obtain  the  optimum  W for  a specified  A.  Using  (2-32)  the 
norm  of  the  gradient  may  be  calculated  for  arbitrary  trial  values  of  A and  W. 
Moreover,  the  norm  of  the  gradient  can  be  expressed  in  terms  of  A alone  if 
W is  set  at  its  optimum  point.  It  can  also  now  be  appreciated  that  the  gradient 
is  not  linearly  related  to  A;  so  that  an  analytic  solution  for  the  optimum  A 
is  not  accessible. 

2.2.4.  Implementing  the  Optimization  Procedure 

As  indicated  in  the  previous  sections  we  can  numerically  attempt  to  mini- 
mize the  MSE  of  our  functional  approximation  by  searching  for  a zero  in  the 
norm  of  thegradient  of  the  MSE  with  respect  to  the  Np  free  parameters  in  the 
hyperplane  coefficient  matrix,  A. 

Any  numerical  solution,  of  course,  is  fundamentally  unable  to  resolve  the 
question  of  whether  or  not  an  absolute  minimum  (as  opposed  to  a rel ati ve 
minimum)  has  been  obtained.  From  a practical  viewpoint,  this  is  not  a severe 
limitation  because,  in  the  final  analysis  the  actual  error  magnitude  will  be 
the  guide  to  acceptance  of  the  approximation. 

However,  problems  can  occur  if  one  attempts  to  blindly  employ  multidi- 
mensional gradient  search  techniques  without  appreciating  the  limitations 
imposed  by  computational  resources. 

Under  constraints  of  finite  computer  resources  it  is  important  to  start 
off  the  numerical  optimization  as  close  as  possible  to  the  optimum  solution 
and  to  make  judicious  tradeoffs  between  suboptimum  minimization  efforts  which 
consume  fewer  computer  resources  but  can  make  significant  progress  towards  a 
good  solution. 
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In  the  next  section  of  this  report  a computational  procedure  is  detailed 
for  implementing  the  optimization  procedure. 

Since  the  procedure  is  a computational  one,  and  since  it  has  been  specif- 
ically implemented  by  software  which  itself  Is  one  of  the  deliverables  under 
this  effort,  the  implementation  description  will  be  unfolded  in  parallel  with 
the  software  description. 
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III.  DETAILED  IMPLEMENTATION  DESCRIPTION 

3.1  OVERVIEW 

The  previous  theoretical  discussion  will  be  specifically  tied  to  a 
software  implementation  in  this  section  of  the  report. 

There  are  Np  P-functions,  each  with  a given  number  of  hyperplanes,  and 
Np  + 1 weights  to  be  determined  so  that  the  mean  squared  error  (MSE)  (defined 
in  Section  2.2.2)  of  the  piecewise-1 inear  fit  is  minimized.  Looked  at  from 
this  point  of  view,  then,  there  results  a minimization  problem.  The  method 
used  to  solve  for  the  best  hyperplanes  and  weights  is  to  employ  a gradient 
search.  This  choice  is  made  because  of  the  nonlinear  dependence  of  the  MSE 
upon  the  hyperplanes  of  the  piecewise-1 inear  fit.  (See  Section  2.2.3.) 

Below  will  be  described  the  normal  sequence  of  operations  in  the  soft- 
ware, assuming  that  the  input  is  data  from  the  Mobility  Model  consistin',  iof 
M Sample  points,  each  one  of  which  has  a fixed  number  n,  of  input  variables 
and  a single  output  variable.  Note  that  the  actual  data  tape  (or  disc  file, 
etc.)  may  have  more  than  n input  variables  per  sample  point,  but  what  is  of 
interest  is  the  actual  number  of  input  variables  selected  from  each  sample 
point,  that  are  to  be  a part  of  the  fit. 

Throughout  this  section  it  is  assumed  that  the  following  are  available: 
the  input/output  data  which  is  to  be  fit,  the  Structural  Parameters  deter- 
minining  the  number  of  P-functions,  number  of  hyperplanes  for  each  P-function, 
and  the  minimum  acceptable  activity  level  (see  Section  3.1.5).  Other  Iteration 
Parameters  (see  Section  3.2.3)  which  determine  the  computational  time  and 
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effort  which  are  to  be  devoted  to  the  fit,  are  also  supplied.  Each  of  the 
P- functions  is  automatically  initialized  by  the  software,  and  then,  using 
the  gradient  search  scheme  referred  to  above,  is  optimized  to  minimize  the 
MSE  within  the  limitations  imposed  on  computer  time  and  effort. 

The  optimization  of  the  hyperplane  coefficients  and  the  weights  occurs 
In  three  distinct  parts:  Marginal  Optimization , Sequential  Optimization 
and  Joint  Optimization.  These  are  described  in  Section  3.1.  It  is  assumed 
in  that  section  that  all  three  will  occur,  as  they  do  in  normal  program 
operation.  Alternative  modes  of  operation  in  which  only  some  of  the  modes 
are  used  will  be  discussed  in  Section  3.2.  Section  3.2  also  discusses  a 
mode  in  which  no  optimization  is  performed  at  all.  Finally,  in  Section  3.3, 
are  discussed  the  principal  building  blocks  of  all  three  types  of  optimization. 
These  are  P- function  initialization,  the  optimizing  of  the  weights  (as  a func- 
tion of  the  hyperplane  coefficients  »nd  the  data),  and  the  gradient  search 
scheme. 

3.1.1  Input  and  Standardization  Of  the  Data 

The  details  of  how  to  input  the  data  are  described  in  Section  4.2, 
but  it  should  be  emphasized  that  the  lexicographic  order  of  the  input  varia- 
bles can  be  very  important  to  the  gradient  search  scheme.  The  reason  for 
this  will  be  evident  in  the  description  of  the  Initialization  and  Marginal 
Optimization  of  the  P-functions  (see  Section  3.2.1).  In  any  case,  the  first 
operation  performed  by  the  software  is  to  put  the  data  into  a suitable 
"standard"  form.  This  standardization  consists  of  modifying  the  input 
variables  so  that  each  one  of  them  has  Its  mean  subtracted  off  and  is 
then  divided  by  Its  standard  deviation.  This  results  in  transformed 
input  variables  which  are  dimensionless.  These  dimensionless  variables 
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are  used  to  eliminate  the  problems  of  scaling  and  to  express  the  input 
variables  in  units  of  standard  deviations  from  the  mean  so  that  the  different 
variables  are  commensurate. 

Formally,  the  transformed  input  variables  {x^,  k = l,...,n>  are  related 
to  the  original  input  variables  x^  by  the  formula, 


(3-1) 


where 


(3-2) 


For  ease  of  notation,  the  x/s  will  be  referred  to  as  x*s  in  the  following 
description  of  the  program  and  its  subroutine.  That  is,  we  shall  always 
assume  that  the  standardized  data  is  being  fitted  and  drop  the  notation. 

Because  a numerical  gradient  search  optimization  is  required,  it  is 
necessary  to  initialize  the  optimization  routine  before  any  optimization 
can  be  performed.  In  addition,  for  computational  efficiency  the  optimization 
is  carried  out  at  a coarse  level  first,  and  proceeds  to  greater  refinement 
in  subsequent  optimization  stages.  It  is  therefore  natural  to  combine  our 
discussion  of  initialization  with  this  first  optimization  level,  which  is 
referred  to  as  Marginal  Optimization. 

3.1.2  Initialization  of  the  P-Functions  and  Marginal  Optimization 

Marginal  Optimization  is  the  first  part  of  the  three  part  optimization 
process  mentioned  earlier.  In  it,  the  following  steps  are  taken:  first  a 
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single  P-f unction  is  initialized  (see  Section  3.2.1)  in  subroutine  INITA. 

Then  weights,  and  w'  are  chosen  in  the  subroutine  OPTW  so  that  w-jP^x.)  + w' 
best  fits  the  data  (see  Section  3.2.2).  Then,  in  the  subroutine  OPT,  a gradient 
search  over  the  hyperplane  coefficients  (see  Section  3.2.3)  of  this  first  P- 
function  and  weights  w-j  and  w'  is  performed  in  OPTW.  In  preparation  for  the 
next  P-function,  the  evaluation  of  w-jP-jfxJ  + w*  is  subtracted  from  the  output 
variable  at  each  data  point  and  entered  in  a residual  array. 

The  steps  for  the  second  P-function  are  the  same  as  those  for  the  first, 
except  that  the  quantity  w2?£ (2L)  + w“  1S  fit  to  the  residuals  remaining  at 
the  end  of  the  marginal  optimization  of  P-|(x)»  Initialization,  initial 
computation  of  the  weights  w2  and  w",  optimization  over  the  coefficients  of 
the  hyperplanes  of  P2(x),  w2,  and  w"  and  a final  recalculation  of  w-j,  w2-  and 
w"  are  performed.  Then  the  residuals 

/ - w^Qc*)  - w2P2(xA)  - w",  l = 1 , . . . ,M 

are  calculated  in  preparation  for  the  third  P-function. 

This  process  continues  until  all  Np  P-f unctions  have  been  involved  one 
at  a time.  At  each  step  of  this  marginal  fitting  procedure,  the  residual 
that  has  not  been  fit  at  the  preceding  step  is  fit  during  the  current  step 
to  the  best  extent  possible  using  a single  P-function  present  alone.  More 
details  pertinent  to  Marginal  Optimization  are  found  In  Section  3.2.1.  The 
resulting  P-function  coefficients  and  weights  are  then  used  as  the  starting 
point  of  the  next  finer  level  of  optimization:  Sequential  Optimization. 
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3.1.3  Sequential  Optimization 

In  Sequential  Optimization,  the  strategy  is  similar  to  Marginal  Optimi- 
zation in  that  we  fit  the  residual  functional  variation  not  yet  accounted  for, 
and  the  fit  is  obtained  by  sequentially  optimizing  one  P-function  at  a time. 
However,  in  this  stage  of  optimization  all  of  the  desired  P-functions  are 
present  while  the  fit  is  adjusted  one  P-function  at  a time.  (This  is  to  be 
contrasted  with  Marginal  Optimization  in  which,  at  any  step,  the  function 
being  fit  to  the  residual  consists  of  a single  P-function.)  The  procedure 
is  as  follows. 

We  start  with  the  P-function  coefficients  and  weights  obtained  from 
Marginal  Optimization.  Then  we  may  assume  that  the  P-functions  P-j,...,P^ 
and  weights  w^  +])  are  available.  The  residual  formed  by  elimin- 
ating P,  : . 


y*  - (w2P2(x*)  + •••  +«Np  PKp  (x*)  + wNp+, ) 


(3-3) 


0 0 

is  calculated  for  each  data  point  (x^  ,y  ) for  X,  = 1,...,M. 

Then  P-jfxJ,  w-j  and  w'  are  optimized  in  the  subroutine  OPT  so  that 


(*.)  + w‘ 


best  fits  the  residual.  This  defines  a new  w^  and  replaces  the  old  value 
of  w^  +-j  with  the  constant  weight  w^  +-j  + w‘. 

The  next  step  is  to  evaluate  the  residuals  formed  by  eliminating  Pg: 


/ - (w^^x1)  + WjPjtx*)  + •••  + wNpPNp<xJt)  + wNp+1)  (3-4) 
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for  l ~ 1,...,M  and  optimize  the  choice  of  the  coefficients  of  Pgfx)*  w2 

and\+v 

This  procedure  continues  sequentially  through  the  P-functions,  as  long 
as  none  of  the  five  stopping  criteria  is  met  (computer  time,  RMS  error,  error 
change,  step  size,  gradient  magnitude  - see  Section  3.2.3).  If  no  stopping 
criteria  has  been  met  at  the  last  P-function,  the  Sequential  Optimization 
will  start  again  with  the  first  P-function  and  continue  sequentially  through 
all  the  P-functions.  The  computer  time  limit  will  force  the  process  to  stop 
if  one  of  the  other  four  criteria  doesn't  stop  it  first.  This  brings  us  to 
Joint  Optimization . > 

3.1.4  Joint  Optimization 

This  third  step  of  optimization  proceeds  .from  the  results  of  the  prior 
optimization  steps.  We  may  therefore  assume  that  P-functions  and  weights 
have  been  previously  defined.  The  procedure  is  simpler  to  describe,  since 
all  the  P-functions  and  all  the  weights  are  simulataneously  changed  in  the 
gradient  search.  The  gradient  search  over  all  the  hyperplane  coefficents 
and  weights  occurs  in  the  subroutine  OPT.  Following  that,  the  subroutine 
OPTW  is  called  to  exactly  determine  the  optimum  weights  for  the  coefficients. 

Joint  optimization  is  a more  ambitious  and  computationally  costly  pro- 
cedure unless  one  is  reasonably  close  to  the  solution  for  the  optimum  hyper- 
plane coefficients.  This  is  because  the  number  of  possible  directions  that 
may  have  to  be  tried  increases  factorially  with  the  number  of  degrees  of 

5 

The  gradient  search  is  computationally  more  efficient  than  repetitive 
numerical  evaluation  of  the  exact  analytical  solution  for  W.  The  analytic 
evaluation  of  W is  therefore  always  reserved  for  the  final  evaluation  of  W 
after  no  further  changes  in  A occur. 
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freedom.  For  this  reason,  the  Marginal  and  Sequential  optimization,  in  which 
the  number  of  degrees  of  freedom  is  much  smaller,  are  used  in  an  attempt  to 
get  a good  starting  point  for  the  Joint  optimization. 

Each  of  these  optimization  stages  has  adjustable  stopping  critera  in 
terms  of  computational  time,  RMS  error,  error  change,  gradient  step  size  and 
norm  of  gradient,  which  are  used  to  control  the  amount  of  effort  spent  in  the 
various  optimization  stages.  (With  the  exception  of  the  time  limits,  which 
must  be  separately  specified  for  each  type  of  optimization,  the  remaining 
stopping  criteria  may  either  be  separately  specified;  or  by  default  the  same 
criteria  is  used  for  all  stages  of  optimization.) 

3.1.5  Removing  Inactive  Hyperplanes 

In  the  course  of  the  optimization,  the  activity  level  of  each  hyperplane 

. * XU 

is  computed.  The  activity  level  of  the  q—  hyperplane  is  the  number  of  sample 
points  on  which  that  hyperplane  is  actually  the  one  on  which  the  maximum  over 
q of  a/*)  • x'  occurs  for  the  i—  P-function.  An  input  supplied  to  the  pro- 
gram is  the  minimum  acceptable  activity  level  per  hyperplane.  If  a hyper- 
plane falls  below  that  level  it  is  eliminated.  This  is  a safeguard  against 
an  overly  ambitious  and  wasteful  specification  of  the  number  of  hyperplanes. 

3.1.6  Returning  to  Original  Coordinates 

At  this  point  a fit  has  been  made  between  the  output  variable  and  the 
standardized  input  variables.  Returning  to  the  notation  of  Section  3.1.1, 

the  Np  P-functions  at  the  end  of  optimization  may  now  be  called  P-|(x_) 

PNp(x)  and  the  weights  may  be  denoted  w-|...,wN  +1  to  represent  the  fact  that 
the  weights  and  P-functions  now  chosen,  operate  on  the  standardized  input 
variables.  The  approximation  achieved  is: 


F(x)  = + ***  + ”nr  ^Np  ® + wNp+l  (3~5) 
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The  program  defines  new  P-functions,  P..(x),  which  operate  on  the  unstandardized 
(i.e.,  original  input)  data  so  that 

Pj(x)  = Pj(£)»  1 = (3-6) 

That  is  the  new  P-functions,  P^,  operating  in  the  un-standardized  input  vari- 
ables give  the  identical  evaluation  as  the  P-functions  coming  from  the  opti- 
mization procedures,  P,.,  operating  on  the  standardized  variables.  It  is 
evident  that  if  we  leave  the  w..  as  they  are  and  substitute  (3-6)  into  (3-5) 
there  results 

t 

Np  Np 

f(£)  - 2 "i  pi& + v+i " pi<^  v+i = F<i) 

i=i  p i=i  p 


Equation  (3-6)  is  attained  by  determining  the  appropriate  new  P-functions 
(i.e.,  hyperplane  coefficients)  which  operate  on  the  unstandardized  input 
variables.  Let  a/^  represent  the  unstandardized  (n  + 1 )-dimensional  vector 
of  hyperplane  coefficients  for  the  q-^  hyperplane  of  the  i—  P-function,  and 
let  denote  the  corresponding  standardized  coefficients  which  resulted 
from  the  optimization  procedures.  The  objective  in  defining  unstandardized 
coefficients  a/^  is  to  have 


+ a 


(D 

q,(n+l) 


+ »(*) 
aq,(n+l) 


(3-7) 
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i.e.,  to  have  the  dot  product  a^  • ^ equal  to  the  dot  product  a/^  • x,*^. 

4 4 

Consider  the  following  relationship  between  standardized  and  unstandardized 

1L  J.L 

coefficients.  For  the  q— hyperplane  of  the  i— P-function  the  unstandardized 
coefficients  are  {a^»  k = l,2,...(n+l)}.  Similar  notation  with  the 
superscript  denotes  the  standardized  coefficients.  Let 


where  and  m^  are  defined  in  Section  3.1,1.  It  can  readily  be  verified 
that  if  the  Equations  of  (3-8)  are  substituted  in  the  left-hand  side  of 
(3-7),  the  right  hand  side  of  (3-7)  results.  Using  (3-8)' the  functional 
approximation  can  be  expressed  in  terms  of  the  unstandardized  hyperplane 
coefficients,  and  they  would  apply  to  unstandardized  independent  variables. 


3.2  ESSENTIAL  INGREDIENTS  OF  THE  THREE  TYPES  OF  OPTIMIZATION 

In  the  previous  section,  the  three  types  of  optimization:  Marginal, 
Sequential,  and  Joint  were  described.  In  Marginal  Optimization,  an  initial- 
ization of  the  P-functions  is  required.  This  ingredient  of  Marginal  Optimi- 
zation is  described  in  Section  3.2.1. 

In  all  three  types  of  optimization,  the  weights  are  optimized  before 
and  after  the  gradient  search  is  performed.  This  ingredient  of  optimization 
is  described  in  Section  3.2.2. 

Finally,  the  gradient  search  itself  is  described  in  Section  3.2.3. 
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3.2.1  P-Function  Initialization  (Used  in  Marginal  Optimization) 

In  Marginal  Optimization  we  attempt  to  fit  the  input/output  data  using 
one  F-function  at  a time  (see  Section  3.1.2).  Ideally,  we  would  begin 
(Initialize)  this  fit  by  specifying  a matrix  of  hyperplane  coefficients  A., 
for  each  of  the  P-functions  {P.(x_);  A.;  i = l,2,...,Np}.  To  specify  the  A. 
we  should  ideally  account  for  the  interdependencies  between  the  n input  vari- 
ables  (the  components  of  £)• 

To  avoid  the  complexities  of  examining  such  interdependencies  we  choose 
an  initialization  procedure  which  can  be  performed  marginally,  i.e.,  one- 

i 

dimension  at  a time.  Thus,  as  we  encounter  each  of  the  P-functions  during 
Marginal  Optimization,  it  will  be  initialized  over  only  one-dimension.  This 
is  equivalent  to  setting  all  but  two  of  the  columns  of  equal  to  zero 

[see  Equation' (2-9)].  One  of  the  two  nonzero  columns  comprises  the  hyperplane 

* . *■  * 

Coefficients  associated  with  the  input  variable  being  initialized  over.  The 
other  nonzero  column  contains  the  required  constant  terms  to  complete  the 
P-function  specification.  This  is  described  more  fully  below. 

The  first  P-function  is  initialized  upon  the  first  input  variable,  the 
second  P-function  upon  the  second,  etc.  If  there  are  more  P-functions  than 
input  variables,  and  if  there  are  n input  variables,  the  (n  + 1)  P-function 
is  also  initialized  on  the  first  input  variable.  In  general,  the  i— P- 
function  is  initialized  on  the  k—  input  variable,  where  i is  equivalent  to 
k mod  n and  k is  between  1 and  n. 

It  Is  Important  to  note  that  this  method  of  initialization  fits  the 

t 

Input  variables  In  the  order  they  are  furnished.  Moreover,  if  there  are 
* 

In  this  case  the  hyperplanes  are  line  segments  and  the  weights  are 
the  slopes  of  the  line  segments. 


fewer  P-functions  than  input  variables  (Np  < n),  then  only  the  first  n input" 
variables  that  are  encountered  will  be  involved  in  the  Marginal  Optimization 
stage.  It  has  been  found  that  better  fits  are  usually  achieved  (particularly 
when  Np  < n)  when  the  independent  variables  are  arranged  in  a particular  order. 
Specifically,  it  has  been  found  that  it  is  useful  to  arrange  the  input  vari- 
ables in  decreasing  order  of  the  absolute  value  of  the  correlation  between 
those  variables  and  the  output  variable.  (See  Section  4.1.1  for  more  explicit 
details.) 

Since  any  single  P-function  presented  alone  is  convex,  and  since  our 
Marginal  Optimization  stage  attempts  the  fit  using  only  one  P-function  at  a 
time,  it  is  reasonable  to  initialize  the  P-functions  so  that  they  start  out 
being  convex.  The  choice  is  somewhat  arbitrary*  and  we  shall  base  our  ini- 
tialization upon  a semicircular  function.  That  is,  our  initialization  will 
be  chosen  so  that  a semicircle  is  approximated  by  a continuous  piecewise- 
1 inear  function  in  which  the  number  of  line  segments  (recall  that  Marginal 
Optimization  initializes  in  one-dimensional  space  where  hyperplanes  are  simply 
line  segments)  of  the  fit  is  equal  to  the  number  of  hyperplanes  selected  for 
the  given  P-function.  The  semicircle  that  is  chosen  has  a radius  of  2 and 
is  situated  as  shown  in  Figure  5.  (The  radius  of  the  semicircle  can  be  kept 
fixed  regardless  of  the  actual  range  of  the  input  variable  because  the  opti- 
mization is  performed  using  the  standardized  form  of  the  input  data 
(Section  3.1.1.)  The  value  of  2 was  chosen  for  the  radius  because  most  of 
the  input  variable  points  should  lie  within  2 standard  deviations  of  the  mean.) 

The  initial  continuous  piecewise-1 inear  approximation  to  the  semicircle 
can  be  described  as  follows.  Suppose  initialization  is  to  be  performed 


Figure  5. 
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Basic  Form  of  Initializing  Function 
(Optimization  Initialized  Over  the  Variable  x^) 


over  the  input  variable  x^.  The  interval  -2  < x^  < + 2 is  divided  into 
NHYP(i)  = Q.  parts,  each  of  length  4/Q.j.  Projecting  these  points,  which  lie 
along  the  x^  axis,  up  to  the  semicircle  defines  Q..  points  on  the  semicircle. 
Connecting  these  points  with  straight  line  segments  provides  the  desired 
initial  fit  to  the  semicircle.  The  resulting  fit  for  Q.  = 4 is  shown  in 
Figure  6.  Mathematically  the  procedure  can  be  described  as  follows: 

Let 

xo  = ~2 

= -2  + 4/Q. 

xq  = -2  + 1 (4/Qt ) 


and  let 

■ -1  + V1  7 *q2  1 

X.  L 

Then  define  the  q-2-  hyperplane's  coefficients  as  follows: 

aqV  “Oifk^iork^n  + l 

aO)  = yq+i  ~yq 

qk  xq+l  “ xq 

a(i)  - Xq+1  yq  ~ yq+l  Xq 

q»"+1  xn+1  - xn 


Figure  6.  Initial  COMPLIAR  Fit  For  Q = 4 


45 


Thus,  for  the  q^- hyperplane,  the  k— coefficient  is  the  slope  between  the 

XL 

q^- and  q + 1 points  on  the  semicircle,  and  the  (n  + 1)  component  is  the  y 
intercept  of  the  line  connecting  these  same  two  points  on  the  semicircle. 

3.2.2  Optimization  of  the  Weights 

Suppose  that  at  a given  point  the  coefficients  for  P-functions  1,...,  Np 
have  been  determined.  The  objective  is  to  determine  the  weights  w-j,  w 2»...» 
wN  +.j  in  (2-13)  so  that  the  MSE  (Equation  (2-15))  is  minimized. 

The  solution  for  the  optimum  weights  is  given  by  Equations  (2-20)  - (2-25). 

3.2.3  Gradient  Search  Scheme 

The  gradient  search  technique  is  used  in  all  three  forms  of  optimization. 

Following  is  a description  of  the  scheme  for  Joint  Optimization.  For  Marginal 

and  Sequential  Optimization,  the  vector,  y (see  (2-17))  consists  only  of  the 

one  P-function's  coefficients  which  are  being  optimized  along  with  two  weights; 

one  is  the  weight  of  the  P-function  being  optimized,  and  the  other  is  a constant 

Weight  (see  Section  3.1.2  and ^3. 1.3).  First  all  the  coefficients  of  all  the 

hyperplanes  of  the  P-function(s)  being  optimized,  along  with  the  weights  are 

arranged  in  one  long  array,  y>  as  in  Equation  (2-17).  Since  the  function 

being  optimized  is  E,  the  gradient  of  E with  respect  to  y is  computed.  This 

-*  — 
gr 

gradient  is  denoted  ~ . Denote  the  initial  set  of  hyperplane  coefficients 

and  weights  y.  Let  Ay  denote  the  change  in  the  value  of  IYqI;  i.e.,  the 

★ 

step  size.  Then  the  trial  value  of  y is 


Y = v m*L  H. 
•Urial  JLo  m 3y 


Remember  that  y is  a vector  in  a space  of  dimension  Np  + Np  + 1. 
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where  m = |3E/3y|.  We  also  define  a step  size  adjustment  factor,  k^  (supplied 
as  an  input  by  the  user).  If  y^ai  actually  achieves  a smaller  error,  E,  then 
it  is  taken  as  the  new  state  variable;  replacing  and  the  procedure  is  re- 
peated. If  the  trial  vector  achieves  a larger  error,  E,  then  the  step  size. 

Ay,  is  divided  by  k^  and  a new  trial  vector  is  calculated.  On  the  other  hand, 
if  a certain  number  (NSMX)  of  decreases  of  E occur  in  a row  (NSMX  is  supplied 
as  an  input  by  the  user),  then  the  step  size  is  multiplied  by  k^.  (The  details 
of  how  the  step  size  is  increased  and  decreased  can  be  found  in  the  IFTRAN 
commented  listing  of  the  subroutine  OPT.  See  Section  6.3.3.,  and  in 
particular.  Figure  11.) 

In  any  case,  the  process  continues  to  iterate  with  five  stopping  criteria 
(listed  in  the  order  they  are  printed  in  the  output)  checked  at  each  iteration: 

(1)  Computer  time 

(2)  RMS  error 

(3)  Error  change  (the  difference  between  the  current 
RMS  error  and  the  previous  RMS  error) 

(4)  Step  size 

(5)  Gradient  magnitude. 

The  stopping  criterion  of  time  must  be  supplied  by  the  user;  the  four  others 
may  be  supplied  and  are  otherwise  given  nominal  values  by  the  program. 

To  complete  the  definition  of  the  gradient  search,  the  calculation  of 
3*E 

the  gradient  -g—  must  be  defined.  The  following  implements  the  procedure 
discussed  In  Equations  (2-27)  - (2-32). 

The  first  Np  elements  of  the  array  y,  are  hyperplane  coefficients,  a ^ . 
The  partial  derivative  of  E with  respect  to  any  a^  may  be  written  (using 
(2-27)  and  (2-31))  as 


V 
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3E 

?T T 

a qk 


M 


Vk 


I E(f(^}  - yA) 
£=1 


if  hyperplane  is 

i.L 

evaluated  on  JL—  point 
otherwise 


The  partial  derivative  with  respect  to  a weight,  w^,  i = 1 Np  is  (using 

(2-20)  through  (2-23)) 

M 

T f = (F(x*)  - /»  P,  CxA)  ; 1 = 1 Np 

1 A=  1 


3E 

3w«i 

Np+» 


2 

M 


M 

E ff  eft 


£=1 


• :■.<  - ' ; ... 

3.3  ALTERNATE  MODES  OF  PROGRAM  OPERATION 

The  normal  sequence  of  program  operation,  where  the  P-f unctions  are 
initialized  and  optimized  by  the  program,  has  been  described  in  Sections  3.1 
and  3.2.  Other  options  of  program  operation  are  available  and  are  described 
below.  They  include: 


• calculate  statistics  on  the  input  data  with  no 
fitting  performed  (3.3.1) 

• refine  previously  determined  P-functions  (3.3.2) 

• calculate  net  hyperplane  statistics  from  previously 
determined  P-functions  (3.3.3). 
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It  is  Important  to  note  that  the  details  of  instrumenting  all  of  the 
alternate  modes  described  in  this  section  are  described  in  Section  4. 

3.3.1  Calculate  Statistics  on  Input  Data 

There  exists  a program  option  to  calculate  the  mean,  standard  deviation, 
maximum,  and  minimum  for  each  of  the  input  variables,  as  well  as  the  output 
variable.  Using  this  same  option,  the  user  also  obtains  the  correlation 
matrix  for  the  collection  of  the  n input  variables  and  the  output  variable. 
This  option  includes  a sub-option  which  provides  a listing  of  all  the  data 
points. 

This  option  was  included  in  the  program  to  allow  the  user  to  do  some 
Initial  checks  on  the  data  and  to  help  the  user  in  the  choice  of  lexicographic 
order  of  input  Variables  (see  Section  3.2.1. on  initialization).  The  normal- 
rule  of  thumb  is  to  order  the  variables  in  decreasing  magnitude  of  their  cor- 
relation  with  the  output  variable.  One  might  wish  to  deviate  from  this  guide- 

M T ‘ * 

\ - 

line,  however,  if,  e.g. , two  input  variables  are  highly  correlated  with  each 
other. 

3.3.2  Refine  Previously  Determined  P-Functions 

Suppose  a given  set  of  P-functlons  is  found  to  give  a reasonably  accept- 
able model  but  it  is  desired  to  use  the  program  to  further  optimize  these 
P-functlons  to  obtain  a better  model.  In  other  words,  suppose  there  exists 
a set  of  P-functions  which  is  to  be  used  as-  a starting  point  In  the  optimi- 
zation scheme.  This  type  of  scenario  may  be  broken  Into  two  distinct  cases. 

In  the  first  case,  the  number  of  P-functions  and  the  number  of  hyper- 
planes In  each  P-functlon  is  kept  the  same,  and  one  simply  initializes  a new 
run  with  these  prior  results.  This  Is  a useful  capability  when  one  desires 
to  examine  the  benefits  of  adding  more  data  points  and/or  allowing  more  time 
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in  the  various  stages  of  optimization.  In  this  first  case,  the  program 
writes  out  the  P-functions  (specifically,  their  hyperplane  coefficients, 
the  weights,  the  number  of  input  variables,  and  the  number  of  hyperplanes  for 
each  P-function)  in  binary  format  at  the  end  of  execution  and  begins  execu- 
tion on  a subsequent  run  using  these  parameters.  The  hyperplane  coefficients 
that  are  available  in  this  first  case  are  in  non-standardized  form. 

In  the  second  case,  the  hyperplane  coefficients  a ^ of  the  P-functions, 
the  weights,  the  number  of  P-functions,  Np,  and  the  number  of  hyperplanes  in 
each  P-function,  Q.,  i=l,...,Np,  are  input  by  the  user  through  the  STPARM 
parameter  list.  The  user  must  furnish  the  hyperplane  coefficients  in  non- 
standardized  form.  Thus,  when  running  from  previously  defined  P-functions 
(in  which  case  Initialization  is  bypassed),  the  program  assumes  that  the  a . 
P-functions  are  input  in  non-standardized  form.  Before  the  program  proceeds 
to  optimize  these  P-functions,  it  puts  them  in  standardized  form.  This  is 
accomplished  using  Equation  (3-8).  The  resulting  P-functions  are  therefore 
standardized  in  the  same  fashion  as  if  automatic  initialization  had  been 
employed. 

An  important  example  of  the  second  case  mentioned  above  occurs  when 
an  attempt  is  made  to  improve  a previously  obtained  fit,  by  trying  an  addi- 
tional P-function.  In  this  situation  one  has  available  the  non-standardized 

* 

P-function  coefficients  of  the  prior  fit  from  the  printed  output  of  the  run. 
These  P-function  coefficients  are  utilized  as  the  initial  coefficients  for 
the  new  fit;  but  the  initial  coefficients  for  the  additional  P-function  must 
also  be  supplied. 

These  are  the  results  contained  in  Table  IB  of  the  output  described 
in  Section  4.4. 
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The  added  P- function  can  be  initialized  any  number  of  ways  by  the  user. 
In  the  absence  of  any  compelling  reason  to  do  otherwise,  the  user  may  choose 
to  manually  Initialize  this  P- function  using  the  same  strategy  as  employed 
in  Automatic  Initialization  (see  Section  3.2.1).  Thus,  suppose  we  wish  to 

i.L 

initialize  the  new  P-function  (call  it  P.)  on  the  k—  variable,  x^.  Suppose, 
also,  that  there  are  to  be  Q..  hyperplanes  In  this  P-function.  Then,  as  in 
Section  3.2.1,  the  initial  values  of  the  coefficients  of  the  Q.  hyperplanes 
comprising  P*  can  be  defined  as  follows  (q=l»...»Qj) 


S qV  * 0.  if  k f 1 or  k f n + 1 


-Ml 

Note  that  these  are  the  desired  Initial  coefficients  , a'^',  of  P^  to  operate 
on  standardized  variables.  The  program,  however,  assumes  that  the  input  hyper- 
plane coefficients  are  not  standardized;  so  we  must  supply  the  coefficients 
of  the  hyperplanes  comprising  P^  so  that,  after  standardization,  they  satisfy 
Equation  (3-9). 


( 
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Consider  setting  the  input  non-standardized  coefficients  a^V  to 

qk 


a^qk  = 0 if  k = 1 or  k ^ n + 1 


a 


(i)  _ 
qk 


aqk/  Qk 


(3-10) 


a 


(i)  , = ~(i)  mk-(1) 

q»(n+l)  aq,(n+l)  ^aqk 


**  / 4 \ . i 

Where  the  values  a'-^7  are  obtained  from  Equation  (3-9)  and  the  and 
are  as  described  in  Section  3.1.6.  Then  it  can  readily  be  verified  that 
after  the  a^  have  been  standardized,  the  result  will  be  the  desired  cocf- 
fi cents  a ^ . 
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IV.  USE  OF  THE  PROGRAM 

This  section  contains  detailed  instructions  on  the  use  of  the  COMPLIAR 
software,  including  a description  of  Inputs  required  and  examples  of  outputs 
furnished  by  the  program. 

Section  4.1  will  describe  how  to  use  the  program  in  a step-by-step 
fashion.  Section  4.2  contains  a complete  list  of  program  inputs  together 
with  their  definitions.  Before  using  the  program,  it  Is  recommended  that 
both  Sections  4.1  and  4.2  be  read  to  get  a complete  idea  of  how  to  use  the 
program. 

Section  4.3  describes  common  problems  which  may  be  encountered  in  program 
operation  and  also  contains  a list  of  error  messages  produced  by  the  program. 
The  printed  output  of  the  program  Is  described  in  Section  4.4.  Finally, 
the  output  oi  P- functions  Is  described  In  Section  4.5. 

4.1  STEP-BY-STEP  PROGRAM  USE 

Below  Is  a description  of  typical  ways  of  using  the  COMPLIAR  software 
to  analyze  data  and  to  produce  the  desired  multidimensional  fit. 

4.1.1  Initial  Data  Analysis  Alone 

By  setting  NPFUNC  = 0,  the  program  will  take  the  input  data  (see  Tables 

4.1  and  4.2)  and  compute  and  print  summary  properties  of  the  data.  The 
minimum,  maximum,  mean,  and  standard  deviation  are  calculated  for  each  of 
the  designated  input  variables  and  the  output  variable.  Also,  the  correla- 
tion matrix  is  calculated  for  the  (n  + l)-dimensional  vector  composed  of 
the  n input  variables  and  the  output  variable.  If  NPFUNC  = -1,  In  addition 
to  the  above,  a listing  of  all  of  the  data  points  is  provided. 
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One  of  the  purposes  of  this  Data  Analysis  option  is  to  familiarize 
the  user  with  the  salient  characteristics  of  the  data.  For  example,  a 
listing  of  the  data  points  permits  identification  of  outliers.  The  cor- 
relation matrix  provides  an  indication  of  the  degree  of  coupling  which 
exists  between  the  various  input  variables  and  also  between  each  of  the 
input  variables  and  the  output  variable.  The  information  contained  in  the 
correlation  matrix  is  important  to  the  working  of  the  program  and  can  in- 
fluence how  good  the  ultimate  fit  will  be  (see  Sections  3.1.2  and  3.1.3). 

f* 

For  example,  for  a given  amount  of  effort  spent  in  the  various  stages  of 
optimization,  it  is  advisable  to  arrange  the  input  variables  in  decreasing 
order  of  importance.  This  is  usually  accomplished  by  arranging  the  input 
variables  in  decreasing  order  of  the  absolute  value  of  their  correlations 
with  the  output  variable.  This  often  works  well,  but  must  be  tempered  by 
one's  knowledge  of  the  problem  at  hand.  It  must’ be  remembered  that  the 
correlation  is  a measure  of  the  linear  association  between  two  variables. 

If  the  correlation  is  near  unity  in  absolute  value,  then  the  two  variables 
are  closely  (and  linearly)  related.  However,  the  correlation  can  be  small 
and  the  variables  can  still  be  closely  (but  not  linearly)  related.  Examples 
can  easily  be  constructed  where  a quadratic  functional  relationship  exists 
between  two  variables,  yet  the  correlation  between  them  is  0.  Therefore, 
if  there  is  even  a strong  belief  that  a variable  is  important,  it  is  useful 
to  place  it  near  the  top  of  the  list. 

Also,  an  examination  of  the  correlation  matrix  may  reveal  a high  cor- 
relation between  two  input  variables,  both  of  which  have  high  correlation 
with  the  output  variable.  In  this  case,  most  of  the  information  should  be 
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obtainable  from  one  of  the  two  variables.  Hence,  a good  strategy  is  to 
include  one  of  them  near  the  top  of  the  list;  but  to  put  the  other  one  at 
the  bottom. 

4.1.2  Initial  Run  of  the  Program 

Having  decided  on  an  order  for  the  Input  variables,  the  next  step  is 
to  run  the  program  and  examine  the  structure  and  quality  of  the  fit  obtained. 
Unless  there  is  a-priori  knowledge  that  the  realtionship  between  input  and 
output  is  convex,  at  least  two  P-functions  should  be  used.  A good  rule  of 
thumb  is  to  try  two  P-functions  with  two  hyperplanes  per  P-function  and 
another  run  with  three  hyperplanes  per  P-function.  After  some  experimenta- 
tion, it  becomes  clearer  how  many  P-functions  and  hyperplanes  per  P-function 
are  necessary  for  a good  fit.  The  activity  vector  may  be  used  as  an  indicator 
of  the  benefit  of  additional  hyperplanes.  If  adding  additional  hyperplanes 
still  leaves  a reasonably  high  level  of  activity  on  most  hyperplanes,  th  n 
the  number  of  hyperplanes  per  P-function  is  probably  appropriate.  However, 
if  many  of  the  hyperplanes  show  up  with  a very  low  activity,  or  even  Q 
activity,  then  the  complexity  of  the  P-functions,  or  their  number,  or  both. 

Is  probably  too  large.  Experimentation  is  required  in  order  to  determine 
reasonable  time  limits  for  the  three  types  of  optimization,  as  well  as 
values  for  the  other  optimization  termination  criteria. 

4.1.3  Starting  Off  the  Optimization  Process  from 
Previously  Determined  P-Functions 

As  outlined  in  Section  3.3.2,  there  are  basically  two  ways  of  starting 
off  the  program  with  previously  determined  P-functions. 

* 

At  least  for  the  amount  of  data  supplied. 
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The  first  is  basically  a s tart/res tart  mechanism  in  the  program:  the 

P-functions,  and  weights  along  with  other  essential  parameters  can  be  saved 
at  the  termination  of  a run  of  the  program  by  specifying  an  appropriate  save 
unit,  ISAVE.  On  another  run,  the  program  can  read  the  old  P- function  from 
unit  IOLD,  where  IOLD  is  identified  as  the  aforementioned  unit.  When  the 
P-functions  and  weights  are  saved  and  restored  in  this  fashion,  they  are  in 
binary  format.  One  may  then  restart  the  program  and  proceed  directly  to  either 
Sequential  Optimization  or  Joint  Optimization.  The  main  utility  of  this  mode 
is  to  enable  one  to  stop  the  program  midway  through  a convergence  to  a P- 
function  and  to  have  the  program  compute  the  hyperplane  statistics.  It  is 
important  to  note  that  when  restarting  the  program  from  P-functions  and 
weights  saved  in  this  fashion,  the  same  P^ function  structure,  i.e.i  number. 
of  P,r functions  and  number  of  hyperplanes  in-each  P-function,  must  be  employed 
as  were  used  in  the  initial  run.  In  fact,  fhe  following  parameters  will  be 

; 1 '■ ' 1 ■ " ''  t ' * ' 

saved  and  restored  when  using  this  option: 

• The  number  of  input  variables 
t The  total  number  of  hyperplanes 
t The  number  of  P-functions 
0 The  number  of  hyperplanes  in  each  P-function 
0 The  hyperplane  coefficients 
0 The  P-function  weights. 

The  second  method  of  starting  from  previously  determined  P-functions 
involves  inputting  the  P-functions  directly  (in  decimal  format)  via  the 
STPARM  parameter  list.  This  option  allows  the  user  to  start  from  a prior 
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fit  and  proceed  with  optimization  while  simultaneously  changing  the  P-function 
structure.  The  details  of  how  to  Initialize  the  new  P- functions  are  given  in 
Section  3.3.2. 

4.2  PROGRAM  INPUT  PARAMETERS 

Below  is  a table  of  the  input  parameters  for  the  COMPLIAR  program. 

They  are  In  the  form  of  parameter  lists.  These  lists  are  described  in  the 

same  order  they  must  be  supplied  to  the  program  and  are  described  in  Tables 

4.1  through  4.5.  The  format  in  each  table  is  the  same:  the  first  column 

contains  the  name  of  the  variable;  the  second  column  contains  a default 

% 

value,  if  one  exists;  or  a "*"  if  a value  must  be  furnished  by  the  user; 
the  third  column  contains  a definition  and/or  explanation  of  the  variable. 


Table  4.1  DD  PARM  Namelist  (Data  Definition) 
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4.3  ERROR  MESSAGES  AND  COMMON  PROBLEMS  ENCOUNTERED 

This  section  describes  some  common  errors  that  can  occur  in  the  use  of 
COMPLIAR  and  techniques  for  fixing  these  errors.  Table  4.6  shows  the  pos- 
sible error  messages,  the  subroutine  from  which  each  is  printed,  and  a brief 
statement  of  how  to  fix  the  error. 

The  error  message  in  OPTW  requires  further  explanation.  It  arises 
if  the  matrix  £*£,  using  the  notation  of  Section  2.2.3,  is  singular. 

The  matrix  £ has  elements  g^.  where,  repeating  Equation  (2-20), 

g£i  = pite*)  » 9&,n+l  = 1 

The  explicit  dependence  of  g^-  on  only  A.,  amongst  all  the  hyperpl ane  coef-„ . 

ficients  is  suppressed  for  the  moment  to  emphasize  the  structure  of  £.  The 

matrix  £*£  is  singular  whenever  the  number  of  linearly  independent  rows  or 

columns  of  £ is  smaller  than  (Np+1).  This  can  occur,  for  example,  if 

M>(Np+l)  (the  usual  case)  and  if  any  two  P- functions  are  identical.  When 

this  happens,  it  is  clear  that  two  columns  of  £ will  be  identical  (linearly 

dependent)  which  makes  £ £ singular. 

In  practice,  there  are  two  common  situations  which  can  cause  two  or 

more  P- functions  to  be  identical  when  OPTW  is  called: 

♦ 

(1)  If  two  or  more  P- functions  are  input  through  the 
STPARM  Parameter  List  with  identical  hyperplane 
coefficients 

(2)  In  the  course  of  P- function  Initialization  and 
Marginal  Optimization  if  too  little  time  is  allotted 
(i.e.,  TIME!  is  too  small)  then  two  P-functions  may 
simultaneously: 


Table  4.6  Error  Messages 
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(a)  have  the  same  values  as  originally  initialized 
(i.e. , no  Marginal  Optimization  actually  takes 
place) 

(b)  have  the  same  number  of  hyperplanes 

(c)  be  initialized  on  the  same  input  variable 
(requires  Np>n). 

The  first  situation  above  can  be  fixed  by  simply  changing  one  of  the  P-functions 
a little. 

The  second  situation  will  arise  whenever  TIME!  is  sufficiently  small  and 
two  P-functions  that  are  initialized  on  the  same  variable  have  the  same  number 
of  hyperplanes.  This  will  happen  when  there  are  more  P-functions  than  input 
variables.  In  this  situation,  if  P.  and  P.  are  two  P-functions,  then  the 

* . v 

initial  valut  of  them  will  be  the  same  if  i = j (mod  n),  where  n is  the 
number  of  ihput  variables,  and  Qy  = Q.  (NHYP(i)  = NHYP(j)). 

* J 

Another  possible  situation  which  can  cause  G^G  to  be  singular  is  if  any 
of  the  P-functions  are  constant.  This  will  cause  two  columns  of  G to  be 
linearly  dependent  because  of  the  column  of  "ones"  present  in  £.  The  fix  is 
obviously  to  remove  these  trivial  P-functions  which  have  zeros  for  the  coef- 
ficients of  the  input  variables.  The  required  constant  will  be  absorbed  by 
one  of  the  remaining  P-functions. 
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4.4  PRINTED  OUTPUT  OF  THE  PROGRAM 

This  section  of  the  report  describes  the  printed  output  produced  by 
COMPLIAR.  The  major  contribution  to  the  output  consists  of  seven  printed 
tables  (referred  to  below  as  Tables  1 through  7 in  accordance  with  the 
example  in  Section  5).  The  reader  will  find  it  helpful  to  read  the  text  of 
this  section  while  glancing  at  the  example  output  contained  in  Section  5 
of  this  report. 

The  printed  output  of  the  program  can  be  readily  divided  into  four 
parts : 

(1)  A listing  of  the  five  parameter  lists  1'n  which  the  user 
specifies  how  the  program  is  to  operate  on  the  data 

(2)  Print  from  each  of  the  phases  of  optimization  used, 
where  the  user  selects  how  much  print  will  be  pro- 

. v'  duced  through  the  absolute  value  of  IPRINT 

(See  Table  4.5) 

1 * * ■ 

(3)  The  series  of  seven  tables  in  which  overall  statistics 
about  the  fit  and  detailed  region  statistics  of  the 
fit  are  displayed  (the  printing  of  these  seven  tables 
occurs  when  | IPRINT | >1 . Table  1 will  be  produced  on 
all  runs,  however. 

(4)  A listing  of  all  the  data  points  used  with  the 
following  information  displayed  for  each  point: 

a)  input  variable 

b)  the  output  variable  * 

c)  the  predicted  output  variable 

d)  the  error 

(The  printing  of  item  (4)  occurs  after  Table  1 and  occurs  when  IPRINT  is 
negative.) 

An  examination  of  the  beginning  of  the  run  in  Section  5 and  Tables 
4.1  through  4.5  shows  clearly  the  first  of  the  four  parts  of  the  output. 


IF 


This  is  the  value  of  the  output  variable  which  is  produced  by  the 
final  repro-model. 
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The  second  phase  of  the  output,  that  which  occurs  during  optimization 
of  the  P- functions,  is  also  exemplified  by  the  portion  of  the  run  in 
Section  5 in  which  optimization  occurs.  The  variables  which  control  the 
amount  of  output  during  optimization  are  IPRINT  and  IPT.  The  elements  of 
the  output  are  decribed  in  Table  4.2.5  in  the  definitions  of  IPRINT  and  IPT. 
These  outputs  are  shown  and  labeled  (i.e. , which  outputs  occur  for  which 
values  of  IPRINT,  since  IPT  is  just  a frequency  designator  for  some  of  the 
print)  at  the  beginning  of  the  P-function  Initialization  and  Optimization 

V 

print  in  the  run  in  Section  5. 

The  main  purpose  of  these  optimization  prints  is  to  give  the  user  some 
indication  of  what  is  actually  going  on  in  the  optimization  process.  It  can 
help  to  answer  questions  like:  are  things  improving  at  a fast  enough  rater 

so  that  on  a subsequent  run  more  time  should  be  allotted?  HowJdoes  adding 
more  P- functions  or  otherwise  changing  the  P-function  structure  affect  the 
quality  of  the  fit  and  the  computational  resources  required?  ' 

The  third  type  of  print,  the  series  of  seven  tables  giving  detailed 
statistics,  is  the  most  important  output  for  analyzing  the  fit  obtained. 

An  examination  of  Tables  1 through  7 in  the  example  of  Section  5 is  recom- 
mended along  with  a reading  of  the  following  description. 

Table  1,  which  is  divided  into  Tables  1A  and  IB,  displays  the  final 
fit  produced  by  the  program:  the  P-functions  appropriate  to  the  unstandard- 

ized (i.e.,  "original"  input)  data.  Table  1 is  printed  after  the  optimi- 
zation is  finished,  hyperplanes  of  insufficient  activity  are  removed 
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(see  Section  3.1.5)  and  the  coefficients  of  the  hyperplanes  are  transformed 

back  to  the  original  (i.e.,  unstandardized)  coordinates.  In  Table  1A,  the 

magnitude  of  the  P- function  weight  has  been  absorbed  into  the  hyperplane 
★ 

coefficients;  whereas  in  Table  IB  the  separate  identity  of  the  hyperplane 
coefficients  and  P-f unction  weights  has  been  preserved.  Table  1A  is  a more 
economical  presentation  of  the  final  fit.  Table  IB  is  included  in  case  the 
user  wishes  to  input  these  P-functions,  or  a modification  of  them  through 
the  STPARM  Parameter  List. 

The  detailed  contents  of  Table  1A  are  as  follows:  FINAL  ERROR  is  the 
value  of  the  overall  RMS  error  obtained  with  the  fit,  and  FINAL  PCTVE  is  the 
corresponding  per-cent  variation  explained. 

The  reader  will  note  the  terminology  "biased"  and  "unbiased"  appended 
to  the  numerical  ERROR  and  PCTVE  results.  This  terminology  is  borrowed  from 
the  general  problem  of  linear  regression  as 'discussed,  for  example,  in  [3]. 
Thus,  if  one  considers  the  expression  for  E in  (2-15)  as  an  estimate  of  the 
variance  of  the  residuals,  E is  a biased  estimate  of  this  variance.  The 
bias  can  be  removed  by  altering  the  constant  factor  used  in  (2-15).  In 
particular,  for  a regression  fit  composed  of  Np  free  parameters,  the  unbiased 
estimate  of  the  variance  of  the  residuals  is  obtained  by  using  the  divisor 
(M-NF)  in  (2-15)  instead  of  M. 

Thus,  when  (2-15)  is  used  as  is,  the  resulting  error  is  referred  to  as 
the  BIASED  MSE,  and  when  the  divisor  (M-Np)  is  used  instead  of  M,  we  refer 

” *Remember that  the  P-function  weights  may  be  either  positive  or  negative. 
It  is  convenient  to  absorb  the  magnitude  of  the  weights  into  the  P-function 
specification  and  separately  keep  track  of  the  algebraic  sign.  The  way  to 
combine  this  information  to  produce  the  final  repro-model  is  illustrated  in 
the  example.  (See  Table  1A  on  page  94.) 

^Draper,  N.  R.  and  H.  Smith,  Applied  Regression  Analysis,  Wiley  & Sons, 
1966,  p.  58-63. 
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to  the  resulting  error  as  the  UNBIASED  MSE.  The  RMS  errors  are  obtained  by 
taking  the  square  root  of  the  appropriate  MSE.  Similarly,  the  appropriate 

* 

value  of  PCTVE  is  obtained  using  the  biased  or  unbiased  estimator  of  the  MSE. 
PCTVE  is  evaluated  from 


PCTVE 


x 100% 


2 2 
where  Op  is  the  biased  or  unbiased  estimate  of  the  MSE  and  a r is  the  variance 

of  the  dependent  (output)  variable  calculated  from  the  data  supplied. 

Several  comments  are  now  in  order  regarding  the  interpretation  of  the 
results.  Generally,  the  biased  RMS  error  and- PCTVE  results  are  satisfactory 
for  judging  the  quality  of  the  fit.  In  any  case,  if  M»Np  there  will  be  very 
little  difference  in  the  numerical  results  anyway  (though  the  biased  RMS  error 
will  always  be  the  smaller  and  thereby  indicate  a'higher  PCTVE).  When  M is 
not  much  larger  than  Np  there  is  a danger  that  too  many  degrees  of  freedom 
are  being  utilized  in  the  fit  for  the  amount  of  data  available.  After  all, 
with  any  finite  set  of  data  points,  if  a sufficient  number  of  parameters  are 
allowed  in  the  fit,  the  data  can  be  fit  exactly.  This  is  "goldmining"  and 
does  not  necessarily  result  in  a good  fit  over  the  overall  input/output 
variable  space.  In  other  words,  the  real  use  of  the  biased  and  unbiased 
numerical  results  is  to  indicate  when  goldmining  or  overfitting  may  be  occurring. 


★ 

Unless  otherwise  noted  in  the  printed  output,  all  RMS  errors  and 
PCTVE  values  are  obtained  using  the  more  common  BIASED  MSE  expression. 
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If  these  numerical  results  are  approximately  equal,  then  M is  sufficiently 
large  for  the  complexity  of  the  fit  attempted.  A large  discrepancy  between 
the  biased  and  unbiased  numerical  results  is  a strong  indication  that  either 
a less  complicated  fit  should  be  used,  or  more  data  must  be  utilized  in  the 
fit. 

As  a further  safeguard,  if  M<Np  (in  which  case  the  unbiased  MSE  would 
produce  a negative  result!)  the  unbiased  RMS  error  output  is  set  to  +999 
and  the  unbiased  PCTVE  is  set  to  -999.  Also,  an  error  message  is  printed 
warning  of  the  paucity  of  data. 

Finally,  we  mention  that  there  are  situations  which  may  cause  PCTVE 
to  be  negative.  This  occurs  if  the  quality  of  the  fit  is  poor.  Instead 
of  printing  the  negative  (and  meaningless)  PCTVE  results,  they  are  redefined 
to  be  zero  to  make  it  easy  to  identify. 

Table  2,  "Overall  Statistics,"  provides  the  mean  and  standard  deviation 
of  all  n Input  variables,  the  output  variable,  and  the  predicted  variable. 
Also  the  full  covaraiance  matrix  and  correlation  matrix  are  provided  for  the 
(n+2)-dimensional  vector  composed  of:  the  n Input  variables;  the  single  out- 
put variable;  and  the  value  of  the  output  variable  predicted  by  the  fit. 

These  matrices  Indicate  whether  or  not  the  correlations  between  each  of  the 
Input  variables  and  the  output  variables  is  similar  to  the  correlation  be- 
tween each  of  the  input  variables  and  the  predicted  variables.  Also,  these 
matrices  reveal  the  correlation  between  the  output  variable  and  the  predicted 
variable. 
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Table  3,  "Net  Hyperplane  Statistics  for  XXX  Regions,"  gives  a listing 
of  regions  in  decreasing  order  of  population.  The  following  information  is 
supplied  for  each  region:  the  region  I.D.,  the  population,  the  percentage 

of  the  total  population,  the  RMS  error  (taken  over  the  points  in  each  region 
separately),  and  the  percent  variation  explained  in  that  region.  While  most 
of  these  terms  are  self  explanatory,  the  region  I.D.  is  not.  This  item  is 
explained  by  the  note  provided  on  the  sample  output  of  Table  3 in  Section  5. 

Table  4,  "Qualitative  Characterization  of  the  Subregions,"  is  designed 
to  give  some  overall  indication  of  how,  in  each  region,  the  mean  of  the 
variable  compares  to  the  overall  mean.  The  unit  of  comparison  is  the 
standard  deviation  of  the  particular  variable.  The  way  the  numbers  in 

i.L 

Table  4 are  arrived  at  is  as  follows:  Call  the  mean  of  the  k~-  variable 

i.L  J.L. 

in  the  j— region  y^  and  let  y^  be  the  mean  of  the  k— variable  over  all 
regions.  Also,  let  represent  the  standard  deviation  of  the  k^-  variable 
overall.  Then  let  t,  '<■<*■ 


kj 


= ^kj  ~ ^k 


The  quantity  R^  is  a measure  of  how  many  standard  deviations  the  mean,  y^ 


/th 


of  the  k — variable  in  region  j,  lies  from  the  global  mean,  y^.  The  value 
printed  in  Table  4 is  derived  from  Rj.  is  as  follows: 
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Value  Printed  for  Region  j and  Variable  k 

-3 
-2 
-1 
1 

2 
i 

3 

Note  that  th?  regions  are  numbered  the  same  In  Tables  3 through  7.  Note  also, 
that  the  variables  on  which  the  calculation  of  Is  done,  include  the  output 
variable  and  the  predicted  variable. 

Table  6,  "Parameters  of  Fit  for  Actual  (Non- Standardized)  Data,"  lists 
the  resultant  hyperplane  coefficients  appropriate  for  the  nonstandardized 
data.  This  table  provides  the  actual  equation,  in  each  subregion,  between 
the  input  and  output  variables. 

Table  7,  "Parameters  of  Fit  for  Standardized  Independent  Variables," 
lists  the  resultant  hyperplane  coefficients  appropriate  to  the  standardized 
Input  variables.  This  display  Is  most  useful  for  guaglng  the  relative  impor- 
tance of  the  Input  variables  in  each  region,  since  these  are  the  coefficients 
of  dimensionless  variables. 


Value  of  R 


W 


RkJ  s -1 


-1  s Rw  s - 1 


-?iRkj<0 
0 -<  Rkj s 1 


ysRkj«1 


1 < R 


kj 
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4.5  OTHER  OUTPUT  OF  THE  PROGRAM 

The  program  also  outputs  P-functions  in  binary  format  when  ISAVE  is 
nonzero,  and  it  outputs  them  to  unit  ISAVE.  As  can  be  readily  seen  from 
inspecting  the  listing  of  CMPLAR,  the  following  are  output  in  the  order 
described 

NDIM1 , NHYPT,  NPFUNC,  NHYP(I),  1=1 NPFUNC 

hyperplane  coefficients 
weights 
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V.  SOME  EXAMPLES 

This  section  of  the  report  is  intended  to  furnish  greater  insight  into 
the  repro-modeling  concept  in  general,  and  the  COMPLIAR  software  in  particular. 

A sample  run  of  COMPLIAR  for  a 4-dimensional  case  is  presented  in  detail 
in  Section  5.1.  This  run  furnishes  a repro-model  which  explains  over  95%  of 
the  variation  between  the  input  and  output  variables. 

Section  5.2  describes  the  results  of  another  run  of  COMPLIAR  for  an 
8-dimensional  case.  The  simplicity  of  the  repro-model,  even  in  this  difficult 
case,  is  a striking  example  of  the  benefits  of  repro-modeling. 

i 

5;i  A SAMPLE  RUN 

Following  is  a sample  COMPLIAR  run.  This  run  utilized  400  data  points 
to  fit  the  dependence  of  vehicle  velocity  (the  ''output"  variable)  upon  four 
terrain  parameters  (the  "input"  variables).  First,  in  Table  5.1,  is  sho  m 
a listing  of  the  input  cards  actually  used  to  generate  this  run.  Following 
this  table  is  the  run  Itself.  The  only  editing  that  has  been  done  in  pre- 
senting the  results  is  to  delete  some  of  the  iteration  print  in  the  opti- 
mization section  to  save  space.  These  deletions  are  Indicated  where  they 


occur. 


Table  5.1  Listing  of  Input  Cards 
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for:  u,$,  army  tank  automotive  research 

AND  DEVELOPMENT  COMMAND  (TARAOCOM) 
. WARREN,  MICHIGAN 
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P-functlon  and  gradient  print  for  iterations  1-7 
have  been  omitted.  
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P-f unction  and  gradient  print  for  iterations  2-15 
have  been  omitted. 
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5.2  AN  8-DIMENSIONAL  EXAMPLE 

The  COMPLIAR  software  was  used  to  furnish  a repro-model  for  an  8- 
dimensional  case.  The  example  utilized  5184  data  points  consisting  of 
values  for  the  8 input  variables  (4  terrain  parameters  and  4 vehicle- 
related  parameters)  together  with  the  associated  values  of  tank  velocity. 
These  data  were  supplied  using  the  '71  Mobility  Model. 

A repro-model  consisting  of  2 P- functions  with  2 hyperplanes  per 
P-functions  wai  capable  of  accounting  for  more  than  94%  of  the  variation 
between  the  Input/output  data. 

The  final  repro-model  is  summarized  in  Table  5.2' below:  Note  the 
simplicity  of  the  result.  The  repro-model  given  in  Table  5.2  can  easily 
be  used  to  furnish  input/output  data  on  a hand  calculator. 
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VI.  PROGRAM  DESCRIPTION 

The  COMPLIAR  program  has  a structure  composed  of  a few  key  elements 
and  several  supporting  elements.  This  structure  becomes  evident  when  a 
top-down  form  of  documentation  is  used.  Top-down  documentation  consists 
of  describing  the  major  functions  of  the  program  and  indenting  at  various 
levels  to  Indicate  the  supporting  structure.  This  type  of  documentation 
is  analogous  to  the  kind  of  outline  used  in  organizing  a text,  where  major 
ideas  and  their  supporting  discussions  are  clearly  evidenced. 

The  program  was  written  In  IFTRAN,  which  is  a TSC  superset  of  FORTRAN. 
IFTRAN  is  a structured  programming  language  and  has  the  facility  to  arrange 
comment  cards  in  the  same  structure  as  the  IFTRAN  program  listing.  When 
appropriate  comment  cards  are  Inserted,  this  produces  what  we  call  a 
"Commented  IFTRAN  Listing."  For  some  subroutines,  the  commented  IFTRAN 
listing  is  an  excellent  way  to  functionally  describe  the  program  structure. 
Where  appropriate,  this  has  been  used. 

Section  6.1  describes  the  main  program.  Section  6.2  describes  the 
routines  used  for  inputting  data  and  other  required  Inputs.  Section  6.3 
describes  the  optimization  routines.  The  statistical  analysis  routines 
are  described  in  Section  6.4.  These  routines  perform  the  analysis  by  regions. 

As  mentioned  earlier,  complete  annotated  listings  of  the  program  can 
be  obtained  from  the  Science  and  Technology  Division  of  the  Tank  Automotive 
Research  and  Development  Laboratory.  These  listings  are  in  FORTRAN  with  the 
comments  of  the  IFTRAN  listing  automatically  Inserted  in  the  appropriate 
positions. 
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6.1  THE  MAIN  PROGRAM  AND  ITS  STRUCTURE 

The  main  program  is  called  CMAIN.  It  allocates  space  for  blank  common 
and  calls  CMPLAR,  which  does  all  the  work  of  the  program.  The  subroutine 
CMPLAR  calls  subroutines  according  to  the  structure  indicated  in  Figure  7. 

The  structure  in  Figure  7 indicates  which  subroutines  are  called  by  the 
various  subroutines.  Figure  7 is  not  a subroutine  tree  as  such.  Figure  8 
shows  the  tree  structure  of  the  program.  Utility  routines  have  been  left 
out  in  both  of  these  figures  in  order  to  simplify  the  description;  a list 
of  them  appear  below.  * 

XMIT  - initializes  arrays 

DOT  - performs  a dot  product 

XNORM  - takes  the  norm  of  a vector 

CODE  - changes  region  I.D.  from  form  used  in  program  ; 
to  more  readable  form  for  output 

TIME  - gets  computer  time* 

The  commented  IFTRAN  listing  of  the  main  subroutine  appears  in  Figure  9. 
It  spells  out  in  detail  how  the  program  works,  from  input,  through  optimiza- 
tion, to  output  of  region  statistics. 

6.2  INPUT  ROUTINES:  SMNPUT  AND  SMINDAT 

The  two  input  routines  are  SMNPUT  and  SMINDAT,  which  are  called  in  that 
order  by  the  main  routine,  CMPLAR.  Their  functions  can  be  readily  determined 
by  examination  of  the  commented  IFTRAN  listings  of  these  subroutines  in 
Figure  10. 

6.3  OPTIMIZATION  ROUTINES 
6.3.1  INITA 

The  initialization  of  the  P-functions  occurs  in  the  subroutine  INITA. 

The  equations  for  this  initialization  are  described  in  Section  3.2.1. 
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I.  Input 

A.  SMNPUT  - processes  4 parameter  lists  (except  DDPARM) 

•1.  DATADF  - processes  DDPARM  parameter  list 
and  defines  INPUT  file 

2.  CORMAT  - computes  correlation  matrix 

II.  Optimization 

A.  Perform  numerical  optimization 

1.  INITA  - initializes  P-functions 

2.  OPTW  - optimizes  weights  for  given 

hyperplane  coefficients 

a.  GAUUSJ  - inverts  G^G  matrix 

3.  OPT  - optimizes  hyperplane  coefficients  and 

weights  using  gradient  search 

a.  Function  subroutine  F - called  every  iteration 

1.  CFCN  - computes  RMS  error,  activity  of 
hyperplanes,  and  gradient  vector 

b.  STOPG  - tests  stopping  criteria 

c.  STPOUT  - output  routine  which  displays 

P-functions  and  other  data 

B.  Output  P-functions 

1.  WSFUNC 

2.  WSFUNF 

C.  Compute  RMS  error  and  percent  variation  explained 

1.  ENERGY  (calls  function  subroutine  F which  calls  CFCN) 

III.  FSTATS  - Analysis  of  Fit  by  Regions 

A.  MNCOV  - computes  overall  statistics 

B.  NETHYP  - manipulates  pointers  and  data  in  order 

to  perform  region  statistics 

1.  HS0RT2  - sorts  arrays 

2.  HS0RT3  - sorts  arrays 

3.  RGSTAT  - computes  region  statistics 

4.  ERCALC  - computes  region  error  statistics 


Figure  7.  Main  Program  and  its  Structure 


CMPLAR 


15 


Figure  8.  Tree  Structure  of  the  Program. 
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SUBROUTINE  CMPUR(A/U,UAA) 

INPUT  DATA  AND  CONTROL  PARAMETERS 
IF  P-PUNCTIONS  ABE  BEING  INPUT  BY  USER 
t IF  P-FUNCTIONS  ARE  IN  BINARY  FORMAT 

, , READ  no,  OF  DIMENSIONS,  NO,  OF  HYPERPLANES,  no.  OP  P-r JNCT10NS, 

, AND  NO,  OF  HYPERPLANES  IN  EACH  P-FUNCTION 

, , IF  INPUT  DIMENSION  DOES  NOT  MATCH  DATA  DIMENSION 

, , , PRINT  ERROR  MESSAGE  AND  STOP 

, . END  IF 

, » RE-ADJUST  POINTERS  TO  BLANK  COMMON  AREAS  FOR  INPUT  P-FUNCTIONS 

, , IF  NOT  ENOUGH  MEMORY  ALLOCATED  FOR  BLANK  COMMON 

, « PRINT  ERROR  MESSAGE  AND  5TOP 

, . END  IF 

, , READ  HYPERPLANES  AND  WEIGHTS 

. END  IF 

, PRINT  INITIAL  P-FUNtTIONS  BEFORE  STANDARDIZATION 

, INPUT  DATA,  STANDARDIZE  DATA 

, INITIALIZE  TEMPORARY  DEPENDENT  VARIABLE  ARRAY  WITH  DEPENDENT 
, VARIABLE  VALUES 
, INITIALIZE  ACTIVITY  AT  HIGH  LEVEL 

, IF  SOME  OPTIMIZATION  IS  TO  BE  DONE  • 

, , standardize  p-functions  • 

. , PRINT  initial  STANDARDIZED  P-FUNCTIONS 

. • , -OPTIMIZE  WEIGHTS 

, , CALCULATE  and  print  RMS  ERROR  AND  PCTVE 

• END  IF 

ELSE  * TIME  HAS  BEEN  ALLOTTED  FOR  MARGINAL  OPTIMIZATION, 

P-FUNCTIONS  TO  BE  INITIALIZED  BY  THE  PROGRAM 
. INPUT  DATA,  STANDARDIZE  DATA 

• SET  UP  TIME  LIMIT  PARAkATERS  AND  PRINT  CURRENT  TIME 

, PRINT  MESSAGE  INDICATING  BEGINNING  OF  STAGE  l OF  OPTIMIZATION 
. INITIALIZE  RESIDUAL  ARRAY  WITH  DEPENDENT  V,  RIABLE  VALUES 

• PRINT  ITERATION  PARAMETER  LIMITS  FOR  STAGE  i 
. DO  FOR  EACH  P-FUNCTION 

, i SET  UP  SIZES  OF  ITERATION  ARRAYS 

. , INITIALIZE  CURRENT  P-FUNCTION  USING  RESIDUAL  ARRAY  FOR  DATA 

, , INITIALIZE  WEIGHTS  FOR  INITIAL  HYPERPLANES  USING  RESIDUAL  ARRAY 

. , AS  DATA 

CALCULATE  RM8  ERROR,  PCTVE  AND  ACTIVITY  . •• 

. . , PRINT  INITIAL  P-FUNCTION  BEING  INITIALIZED 

, . PRINT  INITIAL  ERROR  AND  PCTVE 

, . SET  UP  SIZE  OF  TOTAL  STATE  VECTOR  TO  BE  OPTIMIZED  <2  18  FOR  THE 

, , WEIGHTS) 

, , PERFORM  GRADIENT  SEARCH  OVER  CURRENT  P-FUNCTION  AND  THE  2 WEIGHTS 

. , CALCULATE  AND  PRINT  RMS  ERROR,  PCTVE,  AND  ACTIVITY 

, , PRINT  P-FUNCTION  BEING  fPTIMIZED 

, , OPTIMIZE  WEIGHTS  FOR  EXISTING  P-FUNCTIONS 

, CALCULATE  AND  PRINT  RMS  ERROR  AND  PCTVE 
, , DO  FOR  EVERY  SAMPLE  POINT 

. ...  STORE  RESIOUAL  OF  DEPENDENT  DATA  LESS  THE  EVALUATION  OF 

• . . EXISTING  P-FUNCTIONS  INTO  RESIDUAL  ARRAY 


Figure  9.  Commented  IFTRAN  Listing  of 
Main  Program  (CMPLAR) 


(continued  on  next  page) 
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. END  00  '**' 

END  00 
, RESTORE  VARIABLES 

END  IF  • END  Or  STAGE  l OP  ORtt.MIZATIoN 
IF  SOME  TIME  IS  ALLOWED  FOR  SEQUENTIAL  OPTIMIZATION 
PRINT  ITERATION  PARAMETER  LIMITS  FOR  STAGE  2 

. print  functions  at  start  sequential  optimization  ,, 

SEPEAT  - 

, . , REPEAT  • -> 

, , , set  up  Field  lengths 

, , SAVE  THE  TWO  HyPERPLANE  COEFFICIENTS  WHICH  ARE  TEMPORARILY 

, , OVER-WRITTEN  WITH  2 WEIGHTS 

, , , SET  WEIGHT  OF  CURRENT  P-PUNCTION  BEING  OPTIMIZED  TO  ZERO 

, i , DO  FOR  EACH  POINT  IN  TEMP,  DEp,  VAR,  ARRAV 

t % , , COMPUTE  RESICUAL  OP  DEPENDENT  DATA  LESS  THE  EVALUATION  OF  THE 

,,  , , CURRENT  APPROXIMATING  FUNCTION  (WITH  CURRENTLY  OPTIMIZED  P- 

, , , , FUNCTION  ZEROED  OUT) 

, , END  DO 

, , , SET  UP  INDE'4  LIMITS  FOR  OPTIMIZATION 

, , , INITIALIZE  WEIGHTS  FOR  RESIDUAL  ARRAY 

, , , PERFORM  GRADIENT  SEARCH  FOR  COEFFICIENTS  OF  CURRENT  P-FUNCTION 

, t BEING  OPTIMIZED,  ALONG  WITH  ITS  WEIGHT  AND  CONSTANT  WEIGHT 

, , PRINT  CURRENT  P-FUNCTION 

j.  , SET  WEIGHT  Or  CURRENT  P-FUNCTION  AND  CONSTANT  WEIGHT 

, , , RE8T0RE  VARIABLES  AND  INCREMENT  COUNTERS 

, UNTIL  ALL  P-FUNCTIQNS  HAVE  SEEN  DONE  ONCE  OR  TIME  LIMIT  IS  REACHED 
, RE-INITIALHE  TEMPORARY  DEPENDENT  VARIABLE  ARRAY 
UNTIL  TIME  LIMIT  IS  REACHED,  THEN  PRINT  MESSAGE  INDICATING  ENO  OF 
, SEQUENTIAL  optimization 

, RS-ORTIMIZE  WEIGHTS  FOR  ALL  P-FUNCTI0N8  TOGETHER 

, CALCULATE  AND  PRINT  RMS  ERROR,  PCTVE,  AND  ACTIVITY 
END  IF-  - END  OF  STAGE  2 OF  OPTIMIZATION 
IF  SOME  TIM£  IS  ALLOWED  FOR  JOINT  OPTIMIZATION 
, PRINT  ITERATION  PARAMETER  LIMITS  FOR  STAGE  3 
, PRINT  INITIAL  P-FUNCTIONS 

CALCULATE  AND  PRINT  RM$  ERROR  AND  PCTVE  I • 

, PERFORM  GRADIENT  SEARCH  OVER  ALL  P-FUNCTIONS  AND  WEIGHTS 
, OPTIMIZE  WEIGHTS  FOR  ALL  P-FUNCTIONS  TOGETHER 
, CALCULATE  AND  PRINT  RMS  ERROR,  PCTVE,  AND  ACTIVITY 
, PRINT  FINAL  P-FUNCTJONS  -v' 

END  IF-.  - END  OF  STAGE  $ OF  OPTIMIZATION 
DO  FOR  ALL  HYPERRL ANES  , . -s 

, IF  ACTIVITY  IS  BELOW  MINIMUM 

, , REMOVE  HYPERPLANE 

. END  IP  . ' , f t •:  - 

end  do 

DO  FOR  ALL  INPUT  SAMPLE  POINTS 

, REMOVE  STANDARDIZATION  FROM  INDEPENDENT  VARIABLES  IN  DATA 
END  DO 

IF  HYPERPLANES  ABE  CURRENTLY  STANDARDIZED 
. DO  FOR  iACH  HYPERPLANE 

, , TRANSFORM  HYPERPLANE  COEFFICIENTS  TO  THOSE  FOR  NON-STANOARDIZED, 

, , I,E,  INRUT,  DATA 

, END  DO  - 

END  IP 

IF  SOME  OPTIMIZATION  HAS  BEEN  PERFORMED 
. OPTIMIZE  WEIGHTS  FOR  P-FUNCTION 
END  IF 

CALCULATE  RMS  ERROR,  PCTVE,  AND  ACTIVITY 

PRINT  RMS  ERROR,  PCTVE,  ALL  P-FUNCTIQNS,  AND  ALL  WEIGHTS 

IF  ft  UNIT  HAS  BEEN  SPECIFIED  FOR  SAVING  P-FUNCTIONS 

, SAVE  P-FUNCTIONS  ON  UNIT  ISAVE 

END  IF 

IF  FLAG  INDICATES  NO  FURTHER  ANALYSIS,  STOP 
DO  FOR  ALL  SAMPLE  POINTS 
, COMPUTE  PREDICTED  VALUE  AND  ERROR 

. PRINT  DATA  POINT  • INDEPENDENT  VARIABLES*  DEPENDENT  VARIABLE, 

, PREDICTED  VALUE,  ANO  ERROR 
END  DO 

PERFORM  DETAILED  REGION  STATISTICS 
STOP 

end 


Figure  9.  (Cont.)  Commented  IFTRAN  Listing  of  Main  Program  (CMPLAR) 


SUBROUTINE  SMNPUT 

INITIALIZE  PARAMETERS  AND  ASSIGN  DEFAULTS 
READ  PARAMETERS  IN  FOUR  NAMELISTS 

FILL  OPTIMIZATION  PARAMETERS  with  DEFAULTS  IF  NECESSARY 

READ  DATA  DEFINITION  CARDS 

IF  NPFUNC  IS  NEGATIVE 

, COMPUTE  STATISTICS  FOR  INPUT  DATA 

, STOP 

END  IF 

PRINT  JNRARK,  STPARM,  ITPARM  AND  OTPARM  VARIABLES  WITH  EXPLANATION 

ALLOCATE  BLANK  COMMON 

TEST  SIZE  OF  BLANK  COMMON 

RETURN 

end 


a)  Commented  IFTRAN  Listing  of 
Input  Subroutines  (SMNPUT) 


SUBROUTINE  SHNDAT 
read  data 

FIND  KEAN  AND  STANDARD  DEVIATION  OF  SAMPLES 

NORMALIZE  SAMPLES 

RETURN 

END 


b)  Commented  IFTRAN  Listing  of 
Input  Subroutines  (SMNDAT) 


Figure  10 
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6.3.2  OPTW 

The  optimization  of  the  weights  occurs  in  OPTW.  The  entire  function 

of  this  subroutine  is  to  calculate  the  matrix  £*£  so  that  it  can  be  inverted 

by  GAUSSJ  and  thus  produce  the  optimized  weights.  Since  gTjG  is  symnetric, 

only  the  upper  triangular  part  is  calculated  and  then  the  lower  triangular 

part  is  filled  in,  one  column  at  a time  starting  from  the  first.  The  error 

message  "OPTW  FAILS,  MATRIX  SINGULAR"  is  printed  if  the  subroutine  GAUSSJ 

T 

passes  a value  of  0.0  for  DET,  the  determinant  of  ££.  «■ 

T 

GAUSSJ  calculates  the  inverse  of  the  matrix  £ £ by  using  Gaussian 
elimination.  It  does  not  explicitly  calculate  the  determinant  but  sets  the 
variable  DET  = 0.0  when  it  cannot  find  a pivot  element  which  is  non-zero. 

6.3.3  OPT 

The  subroutine  OPT  performs  a gradient  search  over  the  hyperplane 
coefficients  and  weights  which  are  being  optimized  in  one  of  the  three  types 
of  optimization.  Figure  11  is  a commented  IFTRAN  listing  of  the  subroutine 
OPT  which  explains  how  it  works. 

6.3.4  The  Function  Subroutine  F and  the  Subroutine  CFCN 

The  function  subroutine  F simply  sets  up  index  parameters,  increments 
the  iteration  counter,  Np,  used  in  the  overall  iteration  scheme,  and  calls 
the  subroutine  CFCN  which  gives  the  value  of  the  RMS  error  as  F. 

The  subroutine  CFCN  is  only  called  from  the  function  subroutine  F and 
it  performs  the  following  function: 
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SUBROUTINE  0PT(X,XTE$T,DX,NXX,JJ] 

DESCRIBE  CALLING  PARAMETERS 

SET  UR  ARRAYS  AND  COMMON  AREAS 

OBTAIN  INITIAL  VALUE  OF  ERROR  AND  GRADIENT 

INITIALIZE  ITERATION  COUNTER,  SUCCESSES  COUNTER,  CONSECUTIVE 
SUCCESSES  COUNTER,  AND  ERROR  CHANGE 
LOOP  UNTIL  ONE  OF  THE  ITERATION  TERMINATION  CRITERIA  IS  SATISFIED 
, IF  COMPUTING  TIME  EXCEEDS  ALOTTED  TIME 
, , PRINT  STOPPING  CRITERIA  AND  TIME 

EXIT 

• END  IF 

, IF  STOPPING  CRITERIA  OF  GRADIENT  SIZE,  ERROR  SIZE,  CHANGE  IN 
, ERROR,  OR  STEP  SIZE  IS  MET 

, , PRINT  STOPPING  CRITERIA  AND  TIME 

EXIT 

« END  IF 

, IF  ITERATION  COUNTER  IS  MULTIPLE  OF  PRINT  FREQUENCY  COUNTER 
, PRINT  STATE  VARIABLE  AND  GRADIENT 
, , PRINT  STOPPING  CRITERIA 

, END  IF 

, DO  ELEMENT-WISE  IN  3TATE  VECTOR 

, , DEFINE  TRIAL  STATE  VARIABLE  AS  STATE  VARIABLE  MINUS  THE: 

, PRODUCT  OF  CURRENT  STEP  SIZE  AND  UNIT  GRADIENT 

• END  DO 

, EVALUATE  ERROR  AND  GRADIENT  WITH  THIAL  STATE  VARIABLE 
, IF  ERRUR  WITH  TRIAL  STATE  VARIABLE  IS  GREATER  THAN  ERROR  WIT, 

, CURRENT  STATE  VARIABLE 

, , DIVIDE  CURRENT  STEP  SIZE  BY  FACTOR 

, , RESTORE  GRADIENT  CORRESPONDING  TO  CURRENT  STATE  VARIABLE 

, RE-INITI ALIZE  CONSECUTIVE  SUCCESSES  COUNTER 
, ELSE  TRIAL  STATE  VARIABLE  PRODUCES  SMALLER  ERROR  THAN  CURRENT 
*.  STATE  VARIABLE 

, , SUBSTITUTE  TRIAL  STATE  VARIABLE  INTO  CURRENT  STATE  VARIABLE 

• . S.  SUBSTITUTE  TRIAL  STATE  ERROR  INTO  CURRENT  STATE  ERROR 

J INCREMENT  SUCCESS  AND  CONSECUTIVE  SUCCESSES  COUNTERS 
, , IF  NUMBER  OF  CONSECUTIVE  SUCCESSES  IS  SUFFICIENT 

, , , MULTIPLY  STeP  SIZE  BY  FACTOR 

, , END  IF 

, END  IP 

END  LOOP 

END 


Figure  11.  Commented  IFTRAN  Listing  of  OPT 
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1.  For  each  sample  point  it  evaluates  each  P-f unction, 
and  in  the  course  of  doing  so,  determines  which 
hyperplane  is  evaluated  for  each  P-f unction 

2.  Determines  the  activity  vector  depending  on  the 
result  of  step  1 

3.  Calculates  the  gradient  of  the  hyperplane  coef- 
ficients and  weights  being  optimized 

6.3.5  STPOUT 

The  subroutine  STPOUT  is  an  output  subroutine.  It  first  outputs  one 
line  of  information  as  shown  in  the  listing.  It  then  outputs  the  P-function 
being  optimized  and  its  gradient.  The  logic  controlling  which  prints  occur 
for  which  values  of  IPRINT  and  the  calling  parameter  JJ  can  best  be  seen 
from  the  I FTP AN  listing. 

6.3.6  STOPS 

When  STOPG  is  called  from  OPT,  it  tests  each  of  the  stopping  criteria 
against  their  limits  and  causes  optimization  to  stop  if  any  of  the  four 
stopping  criteria  RMS  error,  error  change,  step  size,  or  gradient  magnitude 
is  exceeded. 

6.3.7  Output  P-Functions:  WSFUNC  and  WSFUNF 

The  two  subroutines  WSFUNC  and  WSFUNF  output  the  P-functions.  In  WSFUNC, 
the  weights  and  the  hyperplane  coefficients  are  output  separately.  In  WSFUNF, 
the  magnitude  of  the  weights  are  multiplied  into  the  hyperplane  coefficients, 
and  the  algebraic  sign  of  the  weights  is  separately  supplied. 

6.3.8  ENERGY 

The  subroutine  ENERGY  is  called  from  CMPLAR  and  is  used  to  determine 
the  RMS  error  and  percent  variation  explained,  both  biased  and  unbiased. 
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It  does  this  by  calling  the  function  subroutine  F (see  Section  6.3.4) 
and  implementing  the  definitions  of  RMS  error  and  PCTVE. 


6.4  FSTATS 

The  subroutine  FSTATS  calls  MNCOV,  prints  out  results  of  the  correla- 
tion coefficient  and  other  variables  calculated  in  MNCOV,  and  then  calls 
NETHYP. 

6.4.1  MNCOV 

The  subroutine  MNCOV  computes  overall  statistics  over  all  the  data 

points.  The  following  are  calculated  for  all  n Input  variables,  the  output 

* 

variable,  and  the  predicted  variable: 


1.  The  means  and  standard  deviations  of  each  variable 

2.  The  covariance  and  correlation  matrices  for  all  n + 2 
variables  together. 


6.4.2  NETHYP 


The  subroutine  NETHYP  does  the  following: 


1.  For  each  sample  point,  it  computes  the  active  hyperplanes, 
sets  up  various  pointers  and  tables  for  sorts 

2.  It  sorts  a table  of  pointers  to  the  data  itself  and 
another  table  of  counters  into  regions  of  Increasing 
numbers  of  points 

3.  Calculates  statistics  for  each  region  and  stores  them 
on  a scratch  file 

4.  Reads  back  the  contents  of  the  scratch  file  for  various 
output  tables. 
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This  subroutine  can  be  readily  understood  by  examining  the  IFTRAN  listings. 
The  comments  are  intended  to  show  how  the  tables  IWORK  and  INET  are  defined 
as  the  program  progresses. 

6.4.3  Statistical  Subroutines  Called  by  NETHYP:  RGSTAT  and  ERCALC 

The  subroutine  RGSTAT  calculates  the  minimum,  maximum,  mean,  and 
standard  deviation  for  all  n input  variables,  the  output  variable,  and  the 
predicted  variable  for  each  region. 

The  subroutine  ERCALC  calculates  the  RMS  error,  the  maximum  error,  ar.d 
the  percent  variation  explained  in  each  region. 


APPENDIX  A 


A-l 

APPENDIX  A 
SUMMARY  OF  NOTATIONS 

COMPLIAR  (continuous  Multivariate  Piecewise-Linear 
Approximati on/Regression) 

Following  Is  a list  of  notations  used  in  this  document.  When  appropri- 
ate, the  corresponding  symbol  used  in  the  program  is  indicated  in  brackets. 

x.  * (x-j. . .x^. . .xn)  vector  input  variable 

x'  = (xv..Xl.. .xnl)  input  variable  vector  with  1 

1 K n augmented  [AAA(IY)] 

y = scalar  output  variable  [AAA (IF)] 

n * number  of  input  variables  [NDIM] 

M = number  of  data  points  [NSAMP] 

i - index  of  data  points  («,  * 1,2,...,M)" 

z_*  = (x4»y£’)  * one  of  the  Input/output  data  points 

Np  a number  of  P*- functions  in  overall  fit  [NPFUNC] 

1 = Index  of  P-function  (i  = 1,2,... Np) 

= number  of  hyperplanes  in  itn  P-function  [NHYP] 

q = index  of  hyperplane 

H^(x.)  = q**1  hyperplane  of  i^  P-function 

overall  matrix  of  hyperplane  coefficients 
for  the  Np  P-functions  [AAA(IHYP)] 


A-2 


a<  = Q 


ra<i>i 

-i 

• 

C?1 — 

.ft) 

-q 

• 

a I1) 
_Q1_ 

matrix  of  hyperplane  coefficients 
for  the  1^  P-function 


hyperplane  coefficients 

4°  * [sqi>"-aqk)"-aq!U)]  for  the  J*  IWPerpUne 

of  the  1xn  P- function 

| Ay | e step  size  In  Gradient  Search  [TK] 

|AyJ  = Initial  step  size  (supplied  by  user)  [input  as  TK, 

0 stored  as  OTK] 

k * step  size  adjustment  factor  [TKFACT] 


' M 

P1  (2L5 Ai ) ■ i<q<Q.  |Hq^(x)|  * 1th  P-function 


W = [Wj,w2,...wN  +^]  s P-f«>nct1on  weight  vector  [AAA(IW)] 
NP 

F(x*iA,W)  B 2w1P1^-4;-1^  + WN  +1  = General  form  of  COMPLIAR  fit 
ifi  P N 

E = MSE  of  approximation  * [y1  - FCx^jA.W)]^ 

1=1 

[function  subroutine  F] 

Nc  * total  number  of  degrees  of  freedom  for  the  fit 
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